
















































FUTURE PUBLICATION OF 
NAVAL RESEARCH LOGISTICS QUARTERLY 


The Office of Naval Research will discontinue publication of the Naval Research Logistics 
Quarterly after distribution of the December 1982 issue. 

The journal will continue to be published, however, and will appear as a nongovernmen¬ 
tal periodical to be published with the cooperation of the Office of Naval Research. 

The new publisher will be John Wiley and Sons, Inc. The Naval Research Logistics Quar¬ 
terly will become a copyrighted journal in the Wiley-Interscience series, and will be sold on a 
subscription basis by the publisher. The cost for a one year subscription will be $60.00 to indi¬ 
viduals and institutions. Wiley-Interscience publication will be initiated with the March 1983 
issue, Vol. 30, No. I. 

A substantial professional discount will be offered to members of the Operations Research 
Society of America and The Institute for Management Sciences. Members of ORSA or TIMS 
may subscribe to the journal for $20.00 per year. 

If you are currently receiving the NRLQ by paid subscription through the Superintendent 
of Documents, U.S. Government Printing Office, and your subscription is due to expire after 
the December 1982 issue, a refund for the unfulfilled portion of your subscription will be 
issued to you by the Government Printing Office. The new publisher has agreed to honor 
requests at the GPO rate for unfulfilled portions which were subscribed to previously. The 
Superintendent of Documents will notify subscribers if they are due refunds because of such 
termination, and refund checks will be mailed subsequently. Subscribers may then elect to 
extend their subscriptions at the GPO rate by contacting the publisher and enclosing their 
endorsed refund checks in payment for the issues due subsequent to termination. An order 
form to facilitate this extension and/or renewal by the new publisher is enclosed. 

Distribution through officially authorized Navy lists and SupDocs dissemination to govern¬ 
ment repository library centers will be discontinued after the December 1982 issue. 

The Office of Naval Research regrets that it cannot provide support incident to such distri¬ 
bution and that those activities requiring copies for official purposes will be expected to bear the 
necessary costs directly. 

Due to restrictions governing the use of GPO subscription lists, the publisher will be 
unable to follow-up with solicitations mailed with listings provided by the Superintendent of 
Documents; hence, future promotional efforts may only be through professional and trade 
sources. Because of this some current recipients may not receive anticipated additional notices 
from the publisher, and may not maintain the continuity of their collections if they rely upon 
the receipt of subscription solicitations from the publisher. Thus, all present recipients who are 
interested in receiving the journal are encouraged to initiate communications with the publisher 
promptly so that their names will not be purged by default. 







It is expected that current editorial policies will be continued in the commercially pub¬ 
lished version. Professor Herbert Solomon of Stanford University will be the new editor in 
charge and will be responsible for the review of manuscripts in the future. New manuscripts 
submitted for consideration for publication may be sent to him immediately (4 copies) at the 
following address: 


Professor Herbert Solomon, Editor 
Naval Research Logistics Quarterly 
Department of Statistics 
Sequoia Hall 
Stanford University 
Stanford, California 94305 

The Office of Naval Research is pleased with the recognition of the contributions to 
research promulgated through the Naval Research Logistics Quarterly. As ONR’s own scientific 
journal it has reflected Navy interests in diverse areas of mathematics research which we 
believed would enrich future work in logistics and systems analysis, and has well served its 
intended purpose of providing a forum for mathematicians and logistics engineers to interact at 
the highest scientific level. The willingness of the publisher and scientific community to per¬ 
petuate the journal on a self-sustaining independent basis is profoundly reassuring that its value 
has been clearly established. With the forthcoming transition to the private sector, the journal 
will be able to better serve the entire scientific community in the area of logistics research. 
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ON SOME STOCHASTIC INEQUALITIES 
INVOLVING MINIMUM OF RANDOM 
VARIABLES* 


Peter Kubat 

The Graduate School of Management 
The University of Rochester 
Rochester, New York 

ABSTRACT 

Let X t be independent IFR random variables and let Y, be independent ex¬ 
ponential random variables such that £|A,1 - £1 K,l for all i - 1.2. n. 

Then it is well known that £| min Uf,)l > £t min < K,)|. Nevertheless, for 
exponentially distributed K,'s and for a decreasing convex function <t>() it is 
shown that 


£[<K min (*,»] < £(<M min (T,))!. 
!</<» !<■■<» 


1. INTRODUCTION 

The bounds for the mean lives of series and parallel systems under various assumptions 
on component life distribution were extensively treated in the well-known book by Barlow and 
Proschan [3]. 

In particular, if we denote by X^Y,) the life length of the Ah component in the first 
(second) series system, then the following proposition holds: 

THEOREM 1: (Barlow and Proschan [3, p. 122]). Let X,(.Y,) have continuous distribu¬ 
tion F,(Gi ) with the mean n,. Let X\ . X„ (Ki, ... , T,) be independent and F, < G,, 

i - 1, .... n, i.e., Ff is star-shaped with respect to G,lF, < G, iff (l/x)G,~ i £,(*) is increasing 
for x > 0]. Then the mean life of a series system using components with lives X\, ... , X n is 

greater than the corresponding system mean life using components with lives Y . Y„, that 

is, 

(1.1) ElminU,, ..., X n )]> £[min(T,. .... Y„)l 

In this note, it will be shown that for F, having Increasing Failure Rate (IFR) and G, 
being exponential we get 

(1.2) £(«Kmin(/Ti. X„)] < £[<t(min(T,. YJl 

•This work was supported in part by the Office of Naval Research under grant No. CNA SUB N00014-76-C-0001. 
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provided that 4>(-) is decreasing convex. A special case of the inequality (1.2) was already dis¬ 
cussed in [2] where bounds for the discounted mean of the series system life was considered. 

2. BOUNDS FOR THE MEAN OF A FUNCTION OF LIFE OF THE SERIES SYSTEM 

Before we prove the main result the following two lemmata are needed. 

LEMMA 1: Let X t — F h /'- 1,2, be two positive continuous r.v.’s such that F, < F 2 
and E[X x \ - E[X 2 ). Then, 

(2.1) £[4>(Y,)] s£ flfcUj)] 
for any convex and bounded function 4»(‘). 

PROOF: For i - 1,2 we have 

£[*(*,)]-■ <P(x)dF i (x) 

' Jo 1 

- -<J> (*)/;(.*•) + f“ 4>(x)F t (x)dx 

where we denoted tlt(x) — <J>'(x) and Fj(x) — 1 — F,(x). Since ♦(■) is bounded then 
-<b(x)F,(x) L - c is finite. Moreover, since <!>(•) is convex then 0(x) is increasing and thus 
all the conditions of Lemma 6.4 [3, p. 1121 are satisfied and 
J 0 " 0(x)£,(x) dx < J o “ <l>(x)F 2 (x)dx. 

The rest of the proof follows easily. □ 

LEMMA 2: Let Z ~ F z (x) - 1 - e~ kx and X ~ *>(*) - 1 - e'** with E[Z] > E\X\. 

Then there exists a unique point, say x 0 , where the density of X crosses the density of Z and 

the crossing is from above. 

PROOF: E[Z\>E[X] implies n > X. Solving If’*' 0 yields x 0 - 

(/* - A) -1 ln(fi/X) >0. □ 

Now we are ready to prove the main theorem. 

THEOREM 2: Let X, ~ F t IFR and Y, - <7, with G, being exponential and £[*,1 - 
E[Y,] for i — 1, ... , n. Let <!>(•) be decreasing convex bounded function on [0,°o). Then 

(2.2) £[4>(min Y,)] < £[«I>(min Y,)l 

PROOF: Denote by X - min(Y,). Let X — F, then F is IFR and there exists an 

exponential r.v. Z with cumulative distribution function (c.d.f.) G such that £[Yl - ElZ). 
Furthermore, F < G [2, p. 1071. 

By the Lemma 1 we have £l<t»(JT)l < £[4>(Z)1, and by Theorem 1 we have 

E\Z\ > £lKl, where Y - min( Y x . K„). Note that Y is also exponential. If 

£IZ] - E[Y) we are done. 
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If E\Z ] > £IK], then there exists a unique point x 0 where the density of Y crosses the 
density of Z and the crossing is from above. Since 

£14>(Z)] - J" <t>(x)\e~ Kx dx 
and 

£fa>( Y)1 - J fl “ 4>( x)ne~»* dx 

we get 

(2.3) £ [<I» (Z) 1 — £[<!>( K)1 — <J>(x) (Xe _Xjf — fie' 11 *) dx 

which can be written as 

(2.4) J* Q l<f»(jr) - <t(x 0 )) (\e~ Xx — Me -MJt ) dx + 4>(x 0 ) J Q (\e~ Kx — tie~ ,ix ) dx. 

The second term of (2.4) is, of course, equal to zero and the first term written as, 

J" o 0 - (<!>(*) - <J>(jr 0 )l (ke~ Kx — fie~ ltx ) dx 

+ (<P(x) - <P(x„)l (ke~ Kx — fA.e~ MX ) dx 4 0. 

is nonpositive, since for x € [0,jr 0 ), 4>(x) - 4>(x 0 ) > 0 (4>(‘) is decreasing) and 

ke~ Kx -ne~ liX < 0 (by Lemma 2). Similarly, for x € (jr 0 .«>), 4 >(jt) - «l»(jr 0 ) < 0 and 
ke~ Kx - ne~ l ‘ x ^ 0. D 


When the function ♦(•) is exponential we obtain the following corollary: 

COROLLARY 1. Let X, — F ( IFR and Y, - G ( exponential, for /- 1. n. 

E[X>] - £[ K,1 - 1//*,. Then for <J>(x) - e~ ax we get 

2 >, 

£{exp{-a(min A^))l < —— 


Finally, let us remark that Theorem 2 will hold even when F) is Increasing Failure Rate 
Average (IFRA), since this condition together with G, being exponential still implies that 
£ < G 13, p. 1071. 

3. EXAMPLES 

Two simple examples will illustrate the theorem. 

EXAMPLE 1: Consider an n-unit series system, where the unit lifetime X, ~ £,. The F, 
is IFR and E[X,] - I/m,. When the system fails it is replaced by a new system. The cost of 
the new system at the time / - 0 is C. The expected discounted cost of replacement at the 

time of system failure is CE[e~ aX 1, where X — min (AT|. X „) denotes the life time of the 

system. 
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It follows from Corollary 1 that the upper bound for the discounted replacement 


C 


^ 2 > 7 ' 


is 


Suppose now that F, is Gamma a, > 1, for all / - 1. rr. i.e., the probaoility 

distribution function (p.d.f.) of X t is 

, \A,x) ar ' - kx 

f(x) - —pp—j— e ' , x > 0, X, > 0. 

In this case F] is IFR, £[Jf;] * a,/X, and thus n, - X/a ( . Moreover, 

CE\e~ aX ] < C/Ia/A + 1), 
where A - £ (Xj/o;). 

EXAMPLE 2: Let X — min(Y|. X„). All X,'s are independent, have IFR distribu¬ 

tion Fi and E[X] - 1//*,. Assume that 4>(x) - (ax + 0) -1 , a > 0, 0 > 0. Clearly, 4>(x) is 
decreasing convex and bounded on (0,oo). ft follows from Theorem 2 that in this case 

£[<MJr)l- £KaX+/3)-'l < r l — Ae~ Ajf dx 
J o ax+/3 

- — e d f°° - e~’~ — £,(</), 

a J * t a 

where A - J/n,, rf- A0/a and £|(</) denotes the Exponential Integral. The numerical 
values of E x (d) are tabulated in [1]. 


REFERENCES 


[1] Abramowitz, M. and I.A. Stegun, Handbook of Mathematical Functions, (Dover Publishing, 
Inc., New York, 1965). 

[2] Agnihothri, S., U.S. Karmarkar, and P. Kubat, "Stochastic Allocation Rules," Operations 
Research, SO, 545-555 (1982). 

[3] Barlow, E.R. and F. Proschan, Statistical Theory of Reliability and Life Testing: Probability 
Models (Holt, Rinehart and Winston, Inc., New York, i 975). 


NAVAL RESEARCH LOGISTICS QUARTERLY 


VOL. 29, NO. 3, SEPTEMBER 1982 








COMPOUND AVAILABILITY MEASURES 


Laurence A. Baxter 

Department of Applied Mathematics ,nd Statistics 
State University of Sew York at Stony Brook 
Stony Brook, Sew York 

ABSTRACT 

An availability measure is the probability that a two-state system modeled 
by an alternating renewal process is available at one or more points or intervals. 
The concept of availability measures is extended to formulae for the joint pred¬ 
iction of availability and numbers of breakdowns (or repairs) of the system 
during a fixed interval. 


1. INTRODUCTION 

Consider a two-state system, i.e., a machine subject to stochastic failure and repair, and 
suppose that the breakdown/reactivation cycle of the machine can be modeled by means of an 
alternating renewal process <e.g., Cox l4l, Chapter 7). Under this assumption, we can derive a 
variety of useful formulae for predicting the reliability of the system. In particular, we can 
determine the distributions of the N(t), the numbers of failures and repairs in a fixed interval 
[1], and we can obtain availability measures , probabilities of the form p{/(t) * I V / f T\ 
where /(f) - 1(0) if the system is operating (failed) at (and where Tis an index set compris¬ 
ing a (finite) series of points or intervals 12], The dependence between the Sit) and the /(/) 
precludes the simultaneous prediction of numbers of failures (or repairs) and availability using 
existing theory. In this paper, we present a series of formulae which generalize the availability 
measures, enabling us jointly to predict numbers of breakdowns (or repairs) and availability. 
We call these formulae compound availability measures. 

Availability measures are of use in assessing the likelihood that a repairable machine 
modeled by an alternating renewal process is available at specified times. If the probability falls 
below a certain threshold, arrangements for securing an adequate back-up can be made. The 
distributions of the Nit), on the other hand, are of use in deciding on the allocation of 
sufficient repair facilities for a fixed time interval. The formulae presented in this paper enable 
us to consider these probabilities simultaneously. 

In Section 2 we examine the simplest case and consider probabilities of the form />{/(/) - 
1, Nit) - j) and p{!iu) - 1 V u € [t,t + x ], Nit) - j). In Section 3 we consider probabili¬ 
ties of the form p[l(t + x) - 1, Nit) — j\, and in Section 4 we consider probabilities of the 
form p[lit) - 1, Nit + x) - Nit) -yj. Some numerical examples for the Weibull/gamma 
process are given in Section 5. It is first necessary to introduce some notation and briefly to 
review some relevant results. 
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Let F and G denote the distribution functions of the failure times and repair times, 
respectively, and let their Stieltjes convolution be denoted 

Kit) - F * Git) - / o ' Fit - u) dGiu). 

Further, the n-fold recursive Stieltjes convolution of K is written K lH> for n — 1,2,3.and 

we define A <ol (r) - 1(0) if t > (<) 0. Let f g and k be the densities corresponding to F , G 
and K , respectively. The renewal counting functions of the numbers of failures and repairs in 
(0,r) are denoted N\it) and N 2 it) iN 2 it) and N 2 it)), respectively, assuming that there is a 
repair (failure) at time 0. The corresponding renewal functions and renewal densities are 
//,(/)- £{A/)(f)) and h,U) - dH t ii)/dt (# - 1,2,3). 

The distributions of the Njit) are as follows: 



| 1 - Fit) 

n — 0 

(1.1) 

p\N\it) - «) “ j F , *<»-■>(,) - F * K {n) it) 

n Js 1 

(1.2) 

p{N 2 it) - n] - *'"'(/) - K ln+ "it) 

n> 0 


|l- Git) 

n - 0 

(1.3) 

p{N}it) - "I “ | G . - c • K ,n, it) 

n > 1 


(see Barlow and Hunter [1]). The point availability of the two-state system is defined as 
A k it) ■ p[l k it) - 1) where k - 0(1) if there is a failure (repair) at time 0. The interval avai¬ 
lability is defined as R k ix,t) - p{l k iu) - 1 V« € [t.t + x)}. It can be shown that 

(1.4) A,(r) - Fit) + F * H 2 it) 

(1.5) A 0 it)- Git) - G • H 2 it) 

(1.6) Riix,t) - Fit + x) + £ h 2 iu)Fit + x - u) du 

(1.7) R 0 ix.t) - h } iu)Fit + x- u)du 
where 

Fit) - 1 - Fit) 121. 

2. COMPOUND MEASURES OF lit) AND Nit) 

The expressions for the simplest compound availability measures, i.e., p{lit) - 1, 
Nit) - j), are of particular interest as surprisingly simple formulae can be obtained for the 
covariances of the /*(/) and the N,it). 


(2.1) 

/>{/.(')- 1 

. N\it) - j\ - 

F • * ( '>(r) 

2 > 0 

(2.2) 

plhit)- 1 

. N 2 it) - j) - 

F • K (j) it) 

j> o 

(2.3) 

p[l 0 it) - 1 

, N 2 it) - j} - 

F • G * K iJ) it) 

j > 0 

(2.4) 

p[l 0 (l) - 1 

. N 3 it) - j) - 

1° 

\ F * G • K ( '- n (f) 

y-o 

J > l. 
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Observe that pllj(r) — 1, bf,(t) = j] - p[l\(t) — 1, N Z U) — y) for all / This is to be 
expected in view of the identity /,(r) — jV 2 (r) - ^(r) + 1. 

The corresponding formulae for P' k (x,t.j) — /?{/*.(«) - 1 V w € [r,r + xl, fy(r) « j) are 
as follows: 


(2.5) 

Pi (x.tj) = • 

P(r + x) 

/„ Fit + x - u) du 

2-0 

2 > 1 

(2.6) 

Pi ( x,t,j) - 

F(t + x) 

k ,J> (u)F(t + x - u) du 

2-0 

j> » 

(2.7) 

Pi (. x.t,j ) - . 

bj(u)F(t + x - u) du 

J > o 

(2.8) 

Pq (x.t.j) - | 

■ 0 

[ J o bj-.\(u) Fit + x - u) du 

2-0 

J > l 


where k {j> (i) - dK {J Hl)/dt and 

J f'gU- u)k { P(u)du j > 1 

V'H.W j — 0 . 


It should be noted that the above formulae, and many of those of Sections 3 and 4, do 
not need to treat the cr / — 0 separately, but this has been done in the interests of clarity. 

The following covariances are readily derived from (2.1M2.4) above. 

(2.9) cov {/,(/).2V,(/)} - A t * H 2 (t) - /4,(r)//,(r) 

(2.10) cov [/,(/), N 2 (i)) - A t • H 2 (l) - A x U)H 2 U) 

(2.11) cov |/ 0 (r),)Vj(/)} - A 0 * H 2 (t) - A 0 U)H 2 (l) 

(2.12) cov {/ 0 (/).Ar 3 (/>) - A 0 * H 2 (t) - A q (i) C// 3 (r) - U. 

These may be obtained directly, in which case the identity 

^jK ,J) U)~H 2 * // 2 (r) + H 2 U) 

will be found helpful, or by exploiting the binary nature of the indicator variable. Thus, for 
example, expression (2.10) follows from the identity 

- 11 - - 01 

and the other covariances may be derived similarly. 
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Observe that the expectation of the product of the indicator variable and N 2 it) is the 
Stieltjes convolution of their expectations. 

3. COMPOUND MEASURES OF lit + x ) AND Nit) 

We now present formulae for probabilities of the form p[l k (t + x) - 1, A'(/) = j). 

(3.1) p{I l (t + x)-\, Ni(t)-j} 

= Fit + x) + f Q fit + .y) A 0 ix - y) dy j = 0 

= So kiJ) (m) ^ ( ' + u)du 

+ /„ / 0 At*-' 1 («)/(/ + y - u) A 0 (x ~ y) dudy 
+ J 0 J* 0 Oy_t (M) git + y — u) ^4 j (at - y) dudy j ^ 1 

where a/D -1 J>t 

l fit ) 7-0 

(3.2) p{/,(f + jr)- 1. 2V Z (/) — y) 

- Fit + x) + J 0 fit + y) <•*<>(* ~ y) dy 

+ J*Jq fiu)git +y - u)A { ix - y) dudy j - 0 
» / 0 k (J) iu)Fit + x - u) du 
+ k {J) iu) fit + y - u) A 0 ix-y) dudy 

+ So So a A u ^ g ^ + y ~ u) A { ix - y) dudy j > 1 

(3.3) p[l 0 it + x) - 1, N 2 it)-j) 

= J Q bjiu) Fit + x - u) dudy 

+ J Q f Q bjiu) fit + y - u) A 0 ix - y) dudy 

+ J* fl k (J) iu) git + y - u) A\ix - y) dudy j > 0 

(3.4) />{/<>(/+ *)- 1. N 3 it)-j) 

= / 0 git + y) Aiix - y) dy 7-0 

- /' bj-iiu) Fit + x- ~ u) du 

+ J 0 bj-\iu) fit + y - u) A 0 ix - y) dudy 

+ Ar 0 ) (u) #(f + y- u) A t ix - y) dudy j > 1. 

Note that, as would be expected, expressions (3.1M3.4) degenerate to (2.1M2.4), 
respectively, at x - 0 and that application of the key renewal theorem shows that 

lim p[l k it + x) - 1, N,it) — 7 } ” —~— p{N/it) — 7 ) 
x— * Ml + M 2 

where ^ and /a 2 are the means of Fand (7, respectively. 
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The corresponding formulae for p[l k iu)~ 1 V u 6 [^+ x,t + x + w], Nft) = j\ are 
obtained from (3.1M3.4) above on replacing Fit + x) by Fit + x + w) and /^(x - y) by 
R k iw,x - y) as appropriate. Further extensions of these formulae to more general compound 
availability measures of the form 

p[l k it\) - / k it 2 ) = • = /*(/„)« 1 , Nj(t) = j) for t < r, < t 2 < • • • < t„ 

are similarly obtained. 

4. COMPOUND MEASURES OF lit) AND Nit.t + x) 

In this section we consider a different type of compound availability measure, that of the 
form p[l k it) - 1, Njit.t + x) - j } where N,it,t + x) * N,it + x) - iIV/it) is the number of 
renewals in [r,/ + xl. 

(4.1) pi/,it)- 1. N/il.t + x) - j) 

= Rxix.t) 7-0 

= J* fl fit + u) p[N 2 ix - u) - j - 1 \du 

+ f 0 f 0 h 2 is) fit - s + u)p(N 2 ix - w) = j - 1 } dsdu j > 1 

(4.2) p[t\(t) - 1, N 2 it,t + x) - j) 

= R x ix,t) + / 0 f(t + u) Gix - u) du 


+ X'X’ h 2 is) fit - s + u) Gix - u) dsdu 7 - 0 

- J Q fit + u)p{N 2 ix - u) - 7 ) du 

+ J Q J o h 2 is) fit - s + u)p[N 2 ix - w) - 7 } dsdu j > 1 

(4.3) p{/ 0 (r) — I, N 2 it,t + x) — /') 

- R 0 ix,t) 7-0 

h 2 is) fit - s + u) p{N 2 ix - u) — 7 - 0} dsdu j > 1 

(4.4) p{/ 0 (/)- 1, A,(/,r + x) - 7 } 

- R 0 ix,t ) + XT h } is) fit - s + u) Gix - u) dsdu 7-0 

h } is)fit - s + m)/»{JV 3 (x - w) - 7 } rfsrfw 7 > 1. 


Observe that, as would be expected, (4.1) and (4.2) degenerate to (1.1) and (1.2), respec¬ 
tively, at t - 0 whereas both (4.3) and (4.4) vanish. Note further that on applying the key 
renewal theorem to the above formulae we see that 

lim p{l\it) - 1, N\it,t + x) - 7 ) - lim p{l 0 it) - 1, N 2 it,t + x) - 7 } 

AVix) 7“ 0 

A J * 0 tyiu) p\N 2 ix - w) - J - 1} du 7^1 
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and that 

lim p[I\(t) - 1, N 2 (l.t + x) - j) - lim p[I 0 (t) - 1, N } (t,t + x) - j\ 

- A f Q \p(u) p[N)(x— u) = j) du j > 0 

where ^4 —= /a.,/(/X] +M 2 >. = F(i)/p ., and ¥(/) = J Q i li(u) du. Thus, we see that the 

asymptotic probabilities are equal to the corresponding probabilities for the equilibrium alternat¬ 
ing renewal process. 

Generalizations to formulae of the form />{/*( u) — 1 V u € lt,t + jc], A/,(/ + w,r + 
w + x) - y'l are readily obtained by obvious modifications to (4.1M4.4). 

5. APPLICATIONS 

One important industrial application of the alternating renewal process is to model the 
sequences of failures and repairs of electricity generating plant, such as boilers and turbo¬ 
generators. The author's experience in analyzing such data shows that, typically, failure times 
are identically distributed, as are repair times, and that these random variables do not appear to 
be dependent. Knowledge of various availability measures is a means of assessing the likeli¬ 
hood that, for example, peak demand will be met over one or more periods. The distributions 
of the renewal counting functions, on the other hand, give some indication of the number of 
breakdowns likely in a fixed interval. These may be used, for example, in deciding the levels 
of stocks of those spare parts which are most commonly used in repair work. The compound 
availability measures enable us to make the two sets of predictions simultaneously. 

As an example, we examine the alternating renewal process with Weibull failure times 
and gamma repair times, the scale parameter being unity in each case, i.e., 

f(t) - at"-' e "*•, git) - ft-* e -'/r(r )). 


Values of 

/>(/,(/) - 1, V,(r) - j) - />{/,(/) - 1. yV 2 (r) - j) - F * K'»(f) 

were evaluated for 0 < t ^ 10 and j - 1,2,3, ... for a variety of values of a and tj. The 
numerical integration was achieved by means of the Cleroux-McConalogue algorithm [3] for 
recursively-defined Stieltjes convolutions using FORTRAN subroutines developed by 
McConalogue [5]. (See McConalogue (6) for further details.) The values of F • 
obtained are illustrated in Figures 1-4 for the following combinations of shape parameters: 

Figure a r\ 

1 1 I 

2 2.5 2.5 

3 4.5 2;5 

4 5 1.25. 
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A COMPARISON OF SEVERAL ESTIMATES OF AVAILABILITY 
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ABSTRACT 

We compare several competing estimates of the availability of a system 
which alternates between two states, "up" and 'down.* in accordance with an al¬ 
ternating renewal process. Both interval and point estimators are compared 
under several special but representative situations The comparison reaffirms 
the validity and robustness of the log-logistic jackknifed estimates However, 
when the point estimates are compared from the intrinsic criterion of probabili¬ 
ty of concentration, the uniformly minimum variance estimate obtained for the 
Markov model performs very well 


1. INTRODUCTION 

In a system which alternates between two stales, "up" or "down," in accordance with an 
alternating renewal process, the long-term point availability (often simply abbreviated as availa¬ 
bility) measures the probability that the process is up at a given instant of time far enough in 
the future. This index measures the long-run expected fraction of lime that the system 
operates satisfactorily. It has been shown in reliability textbooks, (see, for example. Barlow and 
Proschan [l]), that under not too restrictive conditions, this index equals the ratio of the mean 
uptime to the sum of the mean uptime and mean downtime Such two-state models are con¬ 
sidered frequently in the literature in problems related to design and maintenance of computer 
systems, communication systems, power plants, etc System productivity is directly related to 
equipment availability, and. therefore, much attention is currently being placed in the industry 
on the assessment and optimization of plant availability Any availability improvement program 
has to begin with the estimation of current availability In this paper, we examine several con¬ 
tending estimators of availability, and study some of their properties 

In a recent paper, Gaver and Chu |2] have studied the behavior of the jackknife method 
for estimating availability under several different probability models for the underlying distribu¬ 
tions. They have pointed out that frequently in availability studies, it will be difficult to obtain 
clear and unequivocal probability specifications. Their results obtained by Monte Carlo simula¬ 
tion show that the jackknife method is reasonably robust for providing interval estimates and 
valid under several different special but representative probability models. In this paper, we 
carry out an investigation similar in spirit to Gaver and Chu 12). Our point of departure from 
this paper is that we consider an additional estimate of availability (Mazumdar (31), and we also 
consider the intrinsic properties of the estimates from the point of view of point estimation. 
Our main result can be stated as follows: from the point of view of providing interval esti¬ 
mates, the jackknife estimate of Gaver and Chu appears to be suitably valid and robust, and it 
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performs as well or better than the other contending estimates. However, from the point of 
view of point estimation, when the intrinsic property of probability of concentration (Rao [51) 
is considered, the uniformly minimum variance unbiased (UMVU) estimator seems to perform 
as well or better. 

2. THE MODEL AND THE ESTIMATES OF AVAILABILITY 

We assume that the times spent by the process in the "up" and "down" state are mutually 
independent, and we have at hand a sequence of In observations giving the successive up and 
downtimes, (/,. D,, U 2 .D 2 . .... U„, D„ of the system. Except for one case, we shall make 
the usual assumption that the uptimes l/, have the exponential distribution. The downtimes D, 
are identically distributed, and in Section 3, we shall assume that their common distribution 
belongs to the exponential, gamma, lognormal or the Weibull family. The parameter of 
interest, availability, is given by 

(1) A . 

Elt/1 + E[D1 

where E[U ] and £lD] are, respectively, the expectations of U t and D,. 

We consider the following estimators of availability: 

(a) the maximum likelihood estimate: 


(2) 

H 

II 

1 


U + D 

where U 

- £ UJn and D - £ DJn. 

(b) 

the jackknifed maximum likelihood 

(3) 


where 


(3a) 

V - £ Vjn. £ - £ EJn. 

(3b) 

K- nU-(n- 1)Z7_„ 

(3c) 

c --orh) 

(3d) 

E,- nD-Kn- 1 )5_ ( , 

and 


(3e) 

5-, -oHrTT,?, 0 - 

(c) 

the log-logistic jackknifed estimate 

(4) 

^ -_£i- 

*** i+.‘ 
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where 

(4a) Z = - £ Z, 

(4b) Z, - nZ - (n - 1) Z-, 

(4c) Z. ,-ln[I (/,| - In | £ D,j 

(4d) Z - In (17) - In (5). 
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(3) the jackknifed uniformly variance unbiased estimate. 
(6a) A, k . urnvu (/t) - -jj- i ^(/i) 

where 

(6b) W t (n) - M umvll («) — (« — 1) /i- um v U .,(« - D 

A un , vu (n) is given by (Sb), 



(6c) if > 1 

and 

(6d) S_ ( - 


In the terminology related to jackknife estimates, the quantities V,, £„ Z„ W, are referred 
to as pseudovalues. Let 9 be an estimator of the parameter 9 based on a sample of size n. Let 
0, refer to the set of pseudovalues in this situation, / - 1.2. n and let 

(7) 

n 

It has been stated (Miller [4]), that under certain general conditions, the statistic 


( 8 ) 


Vn (9-9) 

-1— T (0, - 9) 2 
n — 1 “f 


has an approximate t distribution with (n - 1) degrees of freedom. Thus, the approximate 
two-sided 100 (I - a)% confidence interval for 9 will be given by (£„, U a ) where 

(9) U a -9 + t a — [ — T (9, - 9) 2 

'-f I "- 1 X 
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In a recent paper, Rao (5] has commented on the wisdom of using minimum mean square 
error as the only criterion for obtaining optimum point estimates. He has suggested that 
because of mathematical convenience and intuitive considerations, the minimum mean squared 
error could be used as a postulate to derive estimators, but he has recommended that their 
acceptability should be judged on more intrinsic criteria such as the probability of concentration 
(PC), where this quantity measures 

(10) PC = Pr (lEstimate — True Valuel < a) 
for different values of a. 

For the purpose of comparing the performance characteristics of the estimates (a)-(e), we 
choose the following three criteria: (i) estimate probability of coverage of the true availability 
using the formula (9) corresponding to a nominal value a; ii) the estimated mean and variance 
of the confidence interval length; and iii) the estimated probability of concentration. (The 
mean squared error was used in [3] to compare the maximum likelihood and the umvu esti¬ 
mate, and we do not repeat the use of this criterion here.) These estimates are obtained by 
performing Monte Carlo simulation experiments where 1,000 synthetic system realizations were 
observed through a cycle of n consecutive up and downtimes. 

Following Gaver and Chu (21, we assume the following distribution forms for the up and 
down times for the purpose of the simulation study. Although these distribution forms 
represent special cases, they are representative of many different practical situations. 

(A) U is exponentially distributed, E\U) - A *; D is exponentially distributed, 

E\D\ - It/,}, {D,| are mutually independent sequences of independent random 

variables. 

(B) U is exponential with E[U) - A -1 ; D is gamma distributed with E[D] - (kfi)'': Var 
\D\ - (y/kfi)~ 2 , where k and n represent, respectively, the shape and scale parame¬ 
ter of the distribution. The successive up and dowtimes are mutually independent as 
in (A). 

(C) U is exponential with E(U) - A" 1 ; D is Weibull distributed with the scale and shape 
parameters y and 8, respectively, given by the following density: 

(11) f D (d) - ; t>0. 

The successive up and downtimes are independently distributed. 

(D) U is exponential with £((/) - A' 1 ; D is lognormally distributed with parameters n 
and <r, given by the following density: 

(12) f D (d ) - - ^1 exp [-1/2(1 wd - M> 2 Ar ? J; d > 0. 
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The successive up and downtimes are independently distributed. 

(E) U is represented by a long-tailed h distribution, see Gaver and Chu (2). For a fixed 
value of the tail-stretch parameter /i, (0 < h < 1), represent an uptime as 

where X is exponentially distributed with E[X] = 1. The sequence {£/j} is one of 
independent random variables, themselves independent of {£>,). The downtimes are 
here assumed independent and exponential. 

In case (A), the ratio U/D is proportional to a statistic having the F distribution with 
(2 n, 2n) degrees of freedom. Thus, an exact confidence interval for the availability parameter 
A can be obtained in this situation with the help of the F-tables. The confidence intervals 
using the F-statistic were also computed. 

3. NUMERICAL COMPARISONS 

In this section we compare the different estimators using the three criteria described 
above for the probability models (A) to (E). 1,000 Monte Carlo realizations were used to 
obtain the estimates of these performance measures. As far as practicable, the same random 
numbers were used for the purpose of the comparison. We examined the case where n - 15. 
Table 1 gives the properties of the confidence intervals obtained from (9). Table 2 gives the 
properties of the point estimates using the PC criterion of equation (10). 

TABLE 1 — Simulation Results' (1000 Runs) Comparing the Two-sided 95% Confidence 
Intervals of Various Estimators of A vailability; n = 15, t = 2.145 
(True Value of Availability ” 0.9901). 


Underlying 

Distributions 

Coverage 

(%) 

Average Length 

Variance of 
Length 

A. (exponential. 

F:9$ 4 

1.65 x 10"* 

3.65 x 10-> 

exponential) 

JK,LL:94.5 

1.87 x 10~ 2 

8.00 x 10- 5 

A - 0.01, fi - 1 

JK.UMVU:92.3 

1.54 x 10“ 2 

4.59 x lO" 5 


JK,MLE;93.2 

1.64 x I0“ 2 

5.02 x I0" 5 

B (exponential. 

F98.4 

1.67 x lO" 2 

2.52 x 10- 5 

gamma) 

JK,LL:94.7 

1.40 x |(T 2 

3.21 x lO" 5 

A - 0.01, M - 3.0 

JK.UMVU93.4 

1.29 x 10- 2 

2.64 x lO' 5 

k - 1/3 

JK.MLE94.8 

1.38 x ir 2 

2.89 x 10" 5 

C. (exponential. 

F:98.0 

1.67 x lO* 2 

2.65 x 10- 5 

Weibull) 

JK.LL:94.3 

1.39 x 10- 2 

3 24 x lO' 5 

A - 0.01. r - M3 

JK,UMVU:92.5 

1.29 x lO" 2 

2.71 x lO" 5 

8 - 2.0 

JK.MLE94.0 

1.38 x 10- 2 

2.98 x ir 5 

D. (exponential. 

F:92.7 

1 68 x lO* 2 

5.30 x 10- 5 

lognormal) 

JK.LL.949 

2.28 x lO” 2 

4.37 x 10-« 

A - 0.01, n - 0.50, 

JK,UMVU:92.5 

1.66 x I0- 2 

8 86 x I0- 5 

<r - 1.0 

JK.MLE898 

1.76 x Hr 2 

9.53 x lO" 5 

E. (long-tailed h. 

F:88 2 

1.77 x I0" 2 

6 74 x I0- 5 

exponential) 

JK,LL:93.4 

2.41 x 10- 2 

1.60 x 10-« 

A - 0.01, M - 1. 

JK.UMVU:92.7 

1.91 x lO" 2 

8.48 x I0" 5 

_ h- 0.2 

JK.MLE:940 

2.01 x lO* 2 

9.18 x lO" 5 


CF: F-stalislic; JK.LL: "Log-Logistic Jackknife;' JK.UMVU "Jackknife UMVU;" 
JK.MLE: "Jackknife Maximum Likelihood Eslimatel 
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TABLE 2 — Simulation Results (1000 Runs) Comparing the Probabilities of 
Concentration for Various Estimates of Availability; n = 15. 

[PC — Pr ({Estimate - Availability < «)] 


r 

Underlying 

E. 

Probability of Concentration (PC) 

,n% 


Distribution 


a - 0.001 

a - 0 0025 

a = 0.005 

a - 01 

A. 

(exponential. 

A* 

229 

538 

83.5 

97,4 


exponential) 

Aft 

20.8 

52.7 

86.1 

98.3 


X - 0.01, n - 1 


22 5 

53.9 

860 

98 3 



Ax*.,. 

21.5 

52.6 

85 8 

98 2 




22.4 

53.4 

834 

97.7 

B 

(exponential. 

A.* 

28 1 

611 

890 

986 


gamma) 

A^ 

266 

61.6 

91.5 

99.0 


X - 0.01, n - 3.0 


28.1 

62.1 

91.7 

99 1 


k - 1/3 


26.9 

61.9 

91.3 

99.0 



Axii 

28.1 

61.8 

898 

98.7 

C 

(exponential. 

Ai* 

262 

59.9 

89.1 

98.3 


Wetbull) 

A /kmlr 1 

25.8 

59.1 

91.8 

98.9 1 


X - 001. y - 1.13 


26.9 

60.8 

924 

99.0 j 


8 - 2.0 

At-... 

26.1 i 

1 59.3 

91.6 

98.9 | 



Axn 

260 

1 59.4 

90.5 

98.4 l 

D 

(exponential. 

An* 

20 5 

1 At i 1 
1 471 1 

82.5 ! 

95 1 


lognormal) 


19.9 

47.4 

83.3 | 

96.6 


X - 0.01. n - 0.50 


19.1 

469 

84.6 

96.4 


<r - 1.0 


20.4 

47.4 

83.6 

96.6 



Axii 

20.0 

47.2 

806 

94.6 

E. 

(long-tailed h. 

Amle 

165 

I 42.6 

73.5 ] 

92.8 


exponential) 

A lk mlt 

15.2 

39.2 

70.8 

94.5 


A - 0.01, n - 1 

A»».» 1 

17.0 

43 4 

74.9 | 

95.1 


8-0.2 1 


155 

400 

71 1 1 

94 1 


J 

■Vii_| 

16.5 

40.7 

72.1 | 

93.0 


Table 1 shows that the log-logistic jackknife estimate provides approximate nominal cover¬ 
age of the true parameter value under a wide variety of situations. In comparison, the jack¬ 
knifed uniformly variance unbiased estimate or the jackknifed maximum likelihood estimate 
does not fare as well. However, when the probability of concentration of the point estimate is 
considered, the uniformly minimum variance unbiased estimate performs very well. In many 
cases, it performs better than the maximum likelihood estimate or the jackknife log-logistic 
estimate. These findings are not surprising in view of the more pronounced locally linear 
characteristic of the log-logistic estimate. 

4. SUMMARY AND CONCLUSIONS 

Several competing estimates of availability of a system which alternates between two 
states — "up” and "down" have been compared. The data is assumed to consist of n complete 
cycles of successive up and downtimes. Monte Carlo simulations carried on several special but 
representative situations indicate the general validity and robustness of the log-logistic jackknife 
estimate of Gaver and Chu. However, when point estimates are considered, from the point of 
view of probability of concentration, the uniformly minimum variance unbiased estimate for 
the Markov case performs very well in these situations. 
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ABSTRACT 

In this paper we consider a simple three-order-statistic asymptotically un¬ 
biased estimator or the Weibull shape parameter r for the case in which all 
three parameters are unknown. Optimal quantiles that minimize the asymptotic 
variance of this estimator, c are determined and shown to depend only on the 
true (unknown) shape parameter value c and in a rather insensitive way. 
Monte Carlo studies further verified that, in practice where the true shape 
parameter c is unknown, using always c with the optimal quantities that 
correspond to c - 2.0 produces estimates, c*. remarkably close to the theoreti¬ 
cal optimal. A second stage estimation procedure, namely recalculating c based 
on the optimal quantiles corresponding to r*. was not worth the additional 
effort. Benchmark simulation comparisons were also made with the best per¬ 
centile estimator of Zanakis (20) and with a new estimator of WyckofT, Bain 
and Engelhard) (181. one that appears to be the best of proposed closed-form 
estimators but uses all sample observations. 

The proposed estimator, c*. should be of interest to practitioners having 
limited resources and to researchers as a starting point for more accurate itera¬ 
tive estimation procedures. Its form is independent of all three Weibull param¬ 
eters and, for not too large sample sizes, it requires the first, last and only one 
other (early) ordered observation. Practical guidelines are provided for choos¬ 
ing the best anticipated estimator of shape for a three-parameter Weibull distri¬ 
bution under different circumstances. 


1. INTRODUCTION 

The cumulative distribution function of a three-parameter Weibull variate A" is given by: 
(1) F(x) - 1 - exp(-[(x - a)/b] c ] x^a 

‘Research performed by this author was supported by the US Office of Naval Research under Contract No. N000I4- 
76-C-0723. 
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where "a" is the threshold, or location parameter, "b " is the scale parameter and V" deter¬ 
mines the shape of the distribution. See Johnson and Kotz, |7| and Mann. Schafer, Sing- 
purwalla [11] for comprehensive information about this distribution and its applications. 

Maximum likelihood estimates (MLEs) of all three parameters have desirable asymptotic 
properties. However, since no closed-form expression for the MLEs exists, the use of a 
numerical iterative procedure is necessary to solve, for a given sample, the corresponding non¬ 
linear optimization problem. This is a computationally difficult problem for which many algo¬ 
rithms fail to provide solutions [19]. Other deficiencies of the maximum-likelihood procedure 
as it applies to this three-parameter distribution are detailed in [16]. 

Practitioners having limited time and access to analysts or computing facilities would thus 
prefer using closed-form (analytic) estimators. However, only rather inefficient estimators cf 
this type for the three-parameter Weibull distribution have been available in the past [11]. 
Recently, more than seventeen closed-form estimators for the three Weibull parameters were 
compared by Zanakis, and the mcst efficient ones identified [20]. 

Of the three Weibull parameters, the shape parameter, c, is particularly important because 
it affects the shape of reliability and hazard rate curves (see [7], [12]), and is the primary cause 
of computational difficulties and estimation errors in MLE problems [10], [21], 

Here, we examine an improved simple percentile estimator, c, of the shape parameter, 
with quantiles determined optimally and independently of the other two Weibull parameters. 

2. A PERCENTILE ESTIMATOR OF SHAPE c WHICH IS INDEPENDENT 
OF THE LOCATION AND SCALE PARAMETERS. 

Let r, 4 t 2 < . •. ^ t„ be the ordered observations from a random sample of size n from 

( 1 ). 


For the rwo parameter Weibull (when the location parameter is known or zero) Dubey [3] 
proposed a two percentile (p, < p k ) estimator for the shape parameter c. 

(2) c = ln(ln(l - p k )/ ln(l - p i ))l\n{i k lt,) 

with i, - t[„ Pi |+i and t k = I np,\ being the greatest integer less than or equal to np,. 

He showed (as was shown earlier by Lieblein [9] in investigating the extreme-value distri¬ 
bution) that its asymptotic variance is minimized when 

(3) p, = 0.16731 and p k = 0.97366. 

Murthy and Swartz [13] derived unbiasing factors and hypothesis tests for c. Recently, Bryson 
[1] argued that the asymptotic efficiency of ? is not especially sensitive to small deviations from 
the optimal percentiles. 

In a recent paper involving simple estimators for all three Weibull parameters [20], five 
estimators of c were compared. These included the following proposed simple percentile esti¬ 
mator, 

(4) c = 0.51n[ln(l - p*)/ln(l - A )]/ln[(r* - r,)/(r, - /,)] 
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where p,, p r p k are three quantiles satisfying 

(5) 0 < p, < p, < p k ^ 1 

(6) —In(1 - pj) = {ln(l - Pj) • ln(l - p k )Y n 

and t,. t k it, < i, < /*) are the corresponding three ordered -,ample observations. (See also 
Johnson and Kotz [7].) 


This estimator requires only three sample observations and is independent of the other 
two Weibull parameters. Thus, when all three parameters are unknown, estimator c is appeal¬ 
ing if the three percentiles can be selected in some optimal way. It was shown before [201 and 
verified here that using (3) and (4), i.e., 

(7) c 2 = 1.494/In [(/* - r,)]/(?, - r,)] 

produces a rather poor behavior. 

The purpose of this paper is to derive optimal (or nearly optimal) quantiles p, and p k ( p, 
then obtained from (6)) for c when all three Weibull parameters are unknown and to compare 
the estimators based on these optimal quantiles with certain othe r s which also require no com¬ 
plicated iterative procedures. Such simple estimators can be used to obtain final estimates or as 
initial estimates for iterative procedures. 


3. OPTIMAL QUANTILES FOR c 


With the use of a theorem of Mosteller [121 and a lemma of Rao [14], the following ana¬ 
lytic expression for the asymptotic variance of cwas obtained: 


>, - Q, 


R, 

1 2 

i , i 

i 2 

or'(Q, - a) 

or' <Qj - a) or' 

Qj ~ Q> 0* - 0 j 

0j-'(ft - 0). 


\Qj-Q> Q>-Qi\\ Or' \ Qj- Q, 


or'or' - qy 2 ) 


R s — pj (1 — p s ) s = i.j.k 


Q , = [-ln(l - Ps )Y u 


so that from (6) 


Qj - (Q t Q„>' /2 - 

Thus, for a given sample size n , this asymptotic variance is a function of the true shape param¬ 
eter c and only two of the three quantiles. 


VOI. 29 , NO 3. SKPTF.MBF.R 1982 


N AVAL RFSFARC H LOGISTIC S QUARTERLY 




422 


S II. /ANAKIS AM) N R. MANN 


A pattern search procedure (61, extended by Zanakis (211 for bounded-variable nonlinear 
optimization problems with or without derivatives, was employed to obtain, for selected values 
of the shape parameter c, the quantiles ft and p k that minimize <6 (ft, p k , c)/n. The results are 
summarized in Table 1 along with the function values at the point given by (3). 


TABLE 1 — Optimal Results for c 


True Shape 
Parameter 

c 

Optimal Values 

Non-Optimal* 

tr<l>(pl,p k ,c) 

Asymptotic 

Inefficiency 

4>{p’,p k ,c)/<f>(p,,p k ,c) 

ft 

P k 

n-<f>(pi,p k .c) 

0.5 

0.0086 

0.9746 

0.230 

0.484 

2.10 

1.0 

0.0048 

0.9816 

1.028 

3.194 

3.11 

1.5 

0.0028 

0.9887 

3.155 

12.314 

3.90 

2.0 

0.0033 

0.9920 

9.096 

34.976 

3.85 

2.5 

0.0051 

0.9932 

23.215 

81.286 

3.50 

3.0 

0.0072 

0.9939 

51.070 

164.374 

3.22 

3.5 

0.0092 

0.9944 

99.545 

300.142 

3.02 

4.0 

0.0109 

0.9947 

176.936 

507.770 

2.87 

4.5 

0.0124 

0.9949 

292.965 

809.210 

2.76 

5.0 

0.0137 

0.9951 

458.762 

1229.437 

2.68 

7.5 

0.0179 

0.9957 

2522.945 

6194.795 

2.46 

10.0 

0.0202 

0.9960 

8314.425 

19586.645 

2.36 


•p'- 0.16731 and pj - 0.97366 


These results reveal that the asymptotically optimal quantiles (p,.p k ) depend on the true 
value of the shape parameter in a rather insensitive way, especially for p k . 

It should be noted that using instead ( p h p k ) from (3), which is optimal for the two- 
parameter Weibull distribution, essentially triples the asymptotic variance of estimator c if all 
three Weibull parameters are unknown. This is more pronounced when the true c - 1.5 to 2.0, 
as the last column of Table 1 indicates. Note also, from Dubey [3] and Kimball (81, that the 
asymptotic variance of c (given by our Equation (2)) when the location parameter, a, is known 
to be or can be set equal to zero, is 0.92 c 2 /n. Thus, incorrectly assuming that "a" is nonzero 
makes the estimation of c much more difficult than necessary. 

4. MONTE CARLO INVESTIGATIONS 

In order to gain further insight into the behavior of the percentile estimator c, two Monte 
Carlo investigations were made. The first investigated the sensitivity of c to small deviations 
from the optimal quantiles and whether a simple iterative scheme could improve the perfor¬ 
mance of c. The second investigation compared c to a very recently proposed efficient percen¬ 
tile estimator, c, based on all order statistics (181 and the best percentile estimator, c\ found in 
our previous study (201. 

Sensitivity of i 

Since the asymptotically optimal quantiles depend upon the true (unknown) shape param¬ 
eter c in sucl¥ an insensitive way, a Monte Carlo study with 1,000 replications was made using 
c, defined by (4), (5), and (6), with ft and p k corresponding to optimal quantile values for 
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c = 2.0 and 3.0. Over ihe range of investigation, i.e., c = 0.5 (0.5)3.5, somewhat better results 
were obtained with c* corresponding to p, = 0.0033, p k = 0.9920 than with one corresponding 
to p, = 0.0072, p k = 0.9939. W e define 

(9) c* - 3.643/ln[(r* - /,)/<», - r,)J 

where /, = /|ooo33*l+i- h ” f|o.<>92n*l+i. and '/ = Onii87*|+i- Note that for most sample sizes of 
interest, r, = r, and t k = i„. 

In this investigation, we also examined the following two estimators: 

c 0 : defined by (4), (5), (6) and based on the optimal combination of p, and p k from 
Table 1, although this would not ordinarily be possible since the true c value is not 
known. For this reason, we also examined the following: 

?**: a single iteration modified estimator 

(10) c** = c(p* p*) 

with p*and pjf selected in a manner prescribed by the value of c*as given in Table 2. 


Table 2 — Rule for Selection of p, and p k in 
Calculation ofc** — c (p*,p*) 


c* 

p* 

Pt 

0.0 - 0.5 

0.0086 

0.9746 

0.5 - 1.0 

0.0048 

0.9816 

1.0-1.5 

0.0028 

0.9887 

1.5 - 2.0 

0.0033 

0.9920 

2.0 - 2.5 

0.0033 

0.9920 

2.5 - 3.0 

0.0033 

0.9920 

3.0 - 3.5 

0.0033 

0.9920 


The schedule shown in Table 2 for selection of combinations of values of p * and p* was 
determined from results involving c 0 and c*. Because of the negative bias of c* for values of 
c > 1.5, a correcting term (ranging from 0.1 to 0.7) was also investigated for calculating c**as 
a function of c*. This tended to decrease the bias of c**al the expense of higher mean square 
errors. Early simulation results also revealed that for n < 100, perfect choice of percentiles 
(those corresponding to the true but unknown c) reduces the bias but not necessarily the MSE, 
particularly when the true c is large. This theoretical estimator, c 0 , has optimal asymptotic pro¬ 
perties and is not necessarily good for small n and large c values. 

Comparison with Two Previous Percentile Estimators c' and c 

As a benchmark comparison, the proposed estimators c* and c** were also compared via a 
similar Monte Carlo study of 5,000 replications with the two best percentile estimators, c' and 
c, suggested in earlier studies [20 and 18]. 

The best of the simple estimators for the shape of a three-parameter Weibull distribution 
examined in [20] was found to be: 
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(11) c'= 2.989/ln[</* - a)/U,- -5)1 
where 

(12) a - (/,f„ - / 2 2 )/(/, + /„ - 2t 2 ) 

was the best estimate of the Weibull location parameter and 
2.989 = In [In (I - 0.97366/ln(l - 0.16731)]. 

Wyckoff, Bain and Engelhardt (18] recently extended their earlier estimator c [4] via the 
following three step procedure: 

First determine the shape parameter initial estimator (C12.20 in 120]) 

H3) c„- 2.989/lnl(/ K , + l -/ l )/(f| )1|>i , +1 - /,)] 

with p h p k given by (3), and the location parameter estimator 

< 14 > a - [r, - 7/» ,/,0 ]/[l - l/ff l f °] 

where t is the sample mean. Note the similarity of estimators (11) and (13). Then the 
Wyckoflf, Bain, Engelhardt estimator is given by: 

(15) c- «*„/[- a) + £ ln(/ f - a)j 

where s - [0.84n] and the constant k„ is tabulated in |4). 

This estimator uses all order statistics and consequently it is not surprising that it was 
found in [18] to have smaller bias and MSE than the much simpler estimator c', which is based 
only on three quantiles. However, c' performed better in [20] than the earlier version of c pro¬ 
posed by Engelhardt and Bain [4]. 

Monte Carlo results comparing the proposed new estimators c* and c** with the previous 
estimators c' and c are shown in Table 3. The nonoptimal estimator c 2 given by (7) was also 
computed in the same simulation, but MSE values were so large and unpredictable for large 
values of c, that no comparisons were necessary. 


TABLE 3 — Comparison of Average Values and Mean Squared Errors ofc',c,c*, and c** 
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Estimator's c m * slightly better performance over c * is not considered worthy of the addi¬ 
tional effort required. Therefore, the latter is instead recommended in all cases—especially 
when the true c is large, but with caution if also n $ 30 

All four estimators behave similarly for small values of <•. particularly as n increases. 
Comparing the proposed estimator c* with our earlier c' we see that the former has a smaller 
bias in almost all cases. The MSE of c* is clearly smaller than that of c’ when n — 100 and 
c > 2.50, In problems with large sample sizes n, c' tends to have almost no variance around 
its expectation; thus, for large n, the MSE of c' consists mostly of the bias squared, with almost 
no contribution from its variance. 

Finally, the bias of c* is about equal to that of c. but the MSE of c is smaller than that of 
c* when c ^ 2.00. Our results about cand c' are consistent with those in (18). 

5. SOME LITERATURE EXAMPLES 

The estimators examined earlier in this paper, along with some other percentile estimators 
used in (20) were applied to seven literature test problems, with known shape parameter. Prac¬ 
titioners may find the results summarized in Table 4 particularly useful. However, these results 
should not be generalized since they represent a static picture to a few examples; something 
like setting your watch at 9:00 and keeping it there: it will be an excellent estimate of time 
twice a day. 




TABLE 4 — Literature Test Problem Results 



Problem 

1 

2 

3 

4 

5 

6 

7 


Source 

15) 

15) 

II7J 

IP) 

12) 

115) 

115) 


Sample Size 

40 

40 

100 

100 

100 

100 

100 



10 

20 



0 

1.975 

0 


5 b 

100 

100 



1 

47.072 

47.072 


H c 

2 

3 

2 

2 

1.2 

1.328 

1.328 


Zanakis c 2 “ C8* 

7.697 

2.027 

1 360 

1.577 

0986 

1.546 

1.344* 

1 

Zanakis <712.214 

2.950 

3.633 

2.001* 

2.139 

1.305 

1.436 

1.600 


Zanakis c'«C12.23t 

1.801 

2.226 

1714 

1.747 

1 080 

1.350 

1.377 


Hassanein C25.234 

1.502 

1.933 

1.755 

1 154 

1 170 

1.321* 

1.041 

5 

Engelhard! 









& Bain C27.23* 

1.657 

1.962 

1.770 

j 1.225 

1.175* 

1.262 

1.139 


Wyckoff, 









Bain & 









Englehardt c 

2.185 

2.571* 

1.970 

1.870* 

1 280 

1.380 

1.444 


Zanakis-Mann c* 

2.026’ 

1.939 

1.712 

1.575 

1 137 

1.231 

1.384 

(p, - 0.0033, p k - 0.9920) 








MLE for c 

2.48 

2.330 

1.783 

1.857 

1.201 + 

1 330+ 

1.767 


•Best percentile estimate for problem. 

+MLE better than best percentile estimate for problem. 
tNotation used in 1201. 


It is interesting to note that in five out of the seven problems, one or more simple per¬ 
centile estimator was more accurate than the MLEs (iteratively obtained on a computer). As it 
was shown earlier [19], [21], this occurs mostly in problems with small shape parameter values, 
particularly when the sample size is small—a fact of particular interest to practitioners having 
limited resources. 


1 


VOL. 29. NO. 3, SEPTEMBER 1982 


NAVAL RESEARCH LOGISTICS QUARTERLY 




426 


S II /AWKIS AMJV K MtW 


6. CONCLUSIONS 

A simple, three-order-statistic estimator c* is developed for estimation of the Weibull 
shape parameter c. when all three parameters are unknown. It was found to compare very well, 
in terms of mean squared error, with an estimator of the same form based on an optimal selec¬ 
tion of order statistics, dependent only upon the unknown true value of c. Compared to our 
previous best simple estimator c' [20], c* has a smaller bias in almost all cases and a smaller 
MSE in problems with large values of n and c. For small values of c. c'and c' are as efficient 
as the new estimator c of Wyckoff, Bain, and Engelhardt [18| which is determined from all 
ordered observations after a single iteration. 

We feel that all three percentile estimators c*, c' and c are useful to practioners and 
researchers, depending on the particular circumstances as suggested in Figure I. 



Figure I. Practical guidelines for choosing an estimator Tor the shape of 
a three-parameter Weibull distribution 


The following comments will further clarify and support our suggestions: 

If a computer code is available, one is naturally tempted to use Maximum Likelihood 
Estimation iterative procedures. However, research has shown that: 
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(i) MLEs sometimes do not exist or represent a local optimum [16]; 

(ii) The shape parameter is the least accurate of the three Weibull parameters to esti¬ 
mate, causing computational difficulties for MLEs as c gets larger [19], 

(iii) In nonregular cases (c < 2), the simple percentile estimator c' was found to be more 
accurate than the corresponding MLE, particularly when n is small [21]; 

(iv) the choice of a percentile estimator like c', c* or c as a starting point may affect the 
speed of convergence, but not so much the final shape parameter estimate of a good 
iterative MLE procedure [21], and 

(v) If the sample size is too large, determination of MLE or even c is very time 
consuming—a real concern in repetitive situations. 

So, even if a computer is available, one may use a percentile estimator initially, addition¬ 
ally or instead. 

It should be noted that practitioners and students (usually nonstatisticians) often need a 
good estimator of the Weibull shape parameter and do not have the time or access to a com¬ 
puter facility, but only a pocket calculator. In such cases, estimator c— which uses all order 
statistics—will be computationally prohibitive for any but extremely small sample sizes. Instead 
we recommend using c* if n ^ 50, c' otherwise. The unusual advantage of c* is that it is 
independent of the other Weibull parameters. 

These quick estimators are also well suited for obtaining point and interval estimates for 
the otherwise unattainable optimum solution to large scale integer, combinatorial or nonconvex 
mathematical programming problems [22]. A sample of heuristic minimum solutions fits natur¬ 
ally a three-parameter Weibull distribution. 

In this paper we also demonstrated how the use of optimization methods can change an 
erratic estimator (c 2 ) into a good estimator (c*), by specifying the appropriate three order 
statistics that minimize the asymptotic variance of this simple estimator. 
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LINEARLY CONSTRAINED PSEUDO-NEWTON METHOD 
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ABSTRACT 

Must Newton-type methods for linearly constrained optimization be either 
of the modified Newton or quasi-Newton variety'’ The contention of this paper 
is that explicity recomputing part of the projected Hessian may be superior to 
both approaches. A computational comparison with MINOS is presented. 


1. INTRODUCTION 

In a previous paper [10), we presented a new algorithm for solving the linearly con¬ 
strained nonlinear programming problem: 

Minimize /(x) 
subject to Ax £ b. 

That algorithm is a modified Newton one, i.e., for a problem involving n decision variables, an 
entire matrix of projected second partial derivatives, which may be as large as n x rt, must be 
computed at each iteration. Such an approach, which requires order n 2 function evaluations 
and order n } arithmetic operations (multiplications and divisions) at each iteration, may not be 
practical when dealing with larger problems, say n > 20. This paper proposes a closely related 
method, which requires substantially less computational work pet iteration, while retaining the 
desirable features of our previous algorithm. 

The chief novelty of our approach in this paper is the idea of explicitly updating part of 
the second derivative matrix at each iteration. To illustrate the idea, the user chooses a "pipe 
width," w, of rows to update each iteration, and the "pipe" moves from upper left to lower right 
from iteration to iteration. For example, if the second derivative matrix is 5 x 5 and 7r - 2, at 
the first iteration the updated elements (indicated by asterisks) would be 

* * *• 0 0 0 

*•000 
0 0 0 0 0 , 

0 0 0 0 0 

0 0 0 0 0 
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at the second iteration, 

0 0**0 
0 0**0 
* * * * 0 

• * * • 0 

0 0 0 0 0 

and at the third iteration 

* 0 0 0 * 

0 0 0 0 * 

0 0 0 0 * 

0 0 0 0 * 


Broadly speaking, Newton-like methods fall into two categories—modified Newton ones 
and quasi-Newton ones. As indicated above, modified Newton methods update an entire 
second derivative matrix at each iteration, while quasi-Newton algorithms iteratively update an 
approximated matrix by adding a matrix of low rank (typically rank 1 or rank 2) to the last such 
approximation. While modified Newton methods generally perform well, they are impractical 
for larger problems. Quasi Newton methods require less work (usually order n function evalua¬ 
tions and order n 2 arithmetic), but are potentially subject to a number of numerical shortcom¬ 
ings, chiefly due to the second derivative approximation process, related to scaling the ele¬ 
ments, restarts, and inability to deal with regions of nonpositive curvature. Our idea is to take 
a middle ground between those two approaches, by letting the user decide how much of the 
second derivative information to recompute each time. Our mathematical analysis shows that 
such an algorithm is a viable idea, and some computational evidence indicates that it may yield 
solutions more rapidly than either equivalent modified or quasi-Newton methods. 

The present method inherits four major strengths from our modified Newton approach. 

First, the theoretical convergence results carry over to the present algorithm. Thus, 
under relatively mild conditions, the algorithm can be shown to converge to a point satisfying 
first order necessary optimality conditions. Under somewhat stronger conditions, the methods 
find a point also satisfying second order necessary conditions. The rate of convergence is also 
shown to be superlinear, that is (||x* +l - x||/l|x* - x||) —* 0 as k — °°, or of order two 
(||x* +1 - x|| < c||x* +l - 3c11 2 for some c), depending on the assumptions made, where (x*) 
is the sequence of points generated by the method and x is the optimal solution. 

Second, the method is reasonably easy to use. It does not require the user to provide 
analytical derivatives. However, while we have had some success in using this approach on 
problems whose derivatives are not smooth (May, Shocker and Sudharshan [11]), the conver¬ 
gence proofs do assume such properties. Our computational experience indicates that the per¬ 
formance of the method is sensitive to only a few parameters, but will converge dependably 
even with highly nonoptimal parameter values. 

Third, the algorithm can deal directly with areas of nonpositive curvature. When minim¬ 
izing a function, and using a Newton-like method, one would like the second derivative matrix 
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(or its projection, if appropriate) to be positive definite, since then the usual search direction 
generation approach yields a descent (improvement) direction. When the second derivative 
matrix is not positive definite, though, a modified approach is required if we want to be sure to 
obtain a useful search direction. Now the second derivative matrix is stored in a factorized 
form, so it can be efficiently used in the solution of the set of linear equations yielding the 
search direction. In our modified Newton method, we utilized the modified LDL T factorization 
of Gill and Murray (6], where L is a unit diagonal lower triangular matrix, and D is a diagonal 
matrix with all positive entries. That factorization, when applied to a nonpositive definite 
matrix A, will actually find the factors of A + £, where E is a diagonal matrix with posiiive 
entries, instead of breaking down, as the usual LDL T factorization would. In our present algo¬ 
rithm, the Gill and Murray LDL T factorization could not be used, since we could not find a 
way of updating it in order n 2 arithmetic when the second derivative matrix is modified by the 
addition of a rank-1 matrix. Instead, we use the factorization of Bunch and Parlett [2], which 
they denote Q T MDM T Q . where Q is a permutation matrix, M unit diagonal lower triangular, 
and D block diagonal, i.e.. it may have entries in the first sub- and super-diagonal, with only I 
x 1 and 2x2 blocks allowed. (We will suppress the permutation matrix Q to simplify the 
notation.) Sorenson [17] has shown how to modify the MDM T factorization when a rank 1 
matrix is added to the matrix under consideration. While modified Newton methods often have 
provisions for dealing with the case of nonpositive definiteness, using special line search tech¬ 
niques, our method appears to be the first one which extends these methodologies to a quasi- 
Newton environment. 

Finally, the algorithm is a "least-constrained" method, in the sense of Lenard [8], which 
means it may find the optimal constraint set, and thus converge, quickly. Since our approach 
uses an e-active constraint strategy, it is desirable to identify those constraints met with equality 
at the optimal solution as quickly as possible. A least constrained method allows for the drop¬ 
ping or adding of several constraints from the active set at each iteration, as opposed to the 
usual simplex-technique related strategy of dropping or adding only one at a time. 


2. THE COORDINATE SYSTEM 


A basic idea in our approach is the representation of the locally feasible region in the 
neighborhood of a point x by a matrix factorization of the active constraint matrix N. That is. 
the columns of N are the normals to the constraint considered active at x. Since we always 
assume that N is full rank, for an n x u matrix N , we can find an n x n orthogonal matrix Q 
and a u x u upper triangular nonsingular matrix R such that 


QU- 


lt 


0 


where 0 is a matrix of zeros. 


Now the last (n — u) rows of @are all orthogonal to the columns of A, which means that 
they define a set of vectors spanning the linear manifold defined by the active constraints. 
Moving away from x along any such row keeps all the active constraint active. Also, the gen¬ 
eralized inverse of N, N + , is given by N + — [/? _l 0]@. The rows of N + have the property that 
movement away from x along the first row of N + , say, corresponds to dropping the first active 
constraint while retaining all the others. 
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The rows of N*. plus the last (n - u) rows of (?, constitute our coordinate system for 
performing an iteration of the method. Since N changes with any change in the active set, so 
does the coordinate system. The method is thus adaptive, adjusting itself to the local structure 
of the feasible region. This adaptability, though, implies a need for proper second derivative 
matrix maintenance, as outlined in Section 5. 

Note that by the orthogonality of Q , the two sets of rows are orthogonal, but that the 
rows of N + are not necessarily mutually orthogonal. 

3. THE ALGORITHM 

The algorithm requires positive real numbers a./3.y.e.r. and p, with p < 1 and 
/3 2 < p/2p. a and /3 are used for move point acceptability, and p is used for a second-order 
Armijo type test for Newton point accep'ability. An e-active constraint strategy is incorporated 
by letting any column A , of A satisfying (A^x — bj < ts be considered active, y is of the 
order of machine precision, and is introduced to avoid numerical difficulties, r is required in 
the theoretical development as a lower bound on the negative curvature curve stepsize but 
could be set arbitrarily small. Also required are an initial feasible point. x n , an initial stepsize 
s° > 0, an initial n-vector <r n of +l’s and -I’s, a symmetric matrix </°, and n, the number of 
Hessian rows / columns to be updated each iteration. 

Let N* and Q be determined by QR factorization of N. Denote the / th column of [A r+ ] r 
by n, and / th column of Q T by </,, that is, 

[A r+ ] r - [w| n 2 ... n u \ and Q T «= [q\ q 2 ... q„\. 

Initialize the current solution point x = x°, the current stepsize s - s°, the current direc¬ 
tion sign indicator tr - o- 0 . the current matrix of QR approximate second partial derivatives 
G — G°, the sequence index k * 0, the last updated column / =■= n , and the local variations 
failure indicator r = 0. 

General Iteration 

Step 1: (Determine active constraints.) Set the square of the gradient norm tj 2 =* 0. 
Determine N for the current x and s. If the active set has changed since the last 
execution of this step, update Q s /V + , and G. 

Step 2: (Computer second order multiplier approximations.) For each constraint /, 

/ - 1,2, .... w, compute A/-, an approximate Karush-Kuhn-Tucker multiplier, by 
evaluating / at two feasible points, a stepsize and a half a stepsize along /?,, and 
using a three point forward difference approximation. If, for any 0, A/, > 0 and 
neither of the two points evaluated along n, yields an improvement, dropping con¬ 
straint i would be undesirable, so set </, - 0. For each / such that A f < 0 or one 
of those two points does yield an improvement, compute an approximate second 
partial derivative A 2 //, set d, - |A/-|, and let tj 2 - r> 2 + (A/)) 2 . 

Let x B denote the point with lowest objective function value amongst all those evaluated 
at the step. 
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Step 3: (First and second derivatives for manifold directions explicitly updated). Let q, 
denote the row of Q corresponding to the j th row of G, which is to be explicitly 
updated this iteration. Then evaluate /at two points, one a feasible stepsize along 
q, away from x and the other a feasible stepsize along —q, away from x. Using a 
central difference formula, compute g t and G ln approximate first and second par¬ 
tial derivatives along q, (the discrepancy in subscripts is due to the fact that we are 
interested in only the last (/» - u) rows of Q ). Set V = « 2 + g 2 . adding g t to the 
norm of the (projected) gradient. 

Update .v B to represent the better of x B from the last step and the lowest function value 
point evaluated at this step. 

Step 4: (First derivatives for manifold directions not updated.) Let q, denote the row of Q 
corresponding to the j th row of G, which is not to be explicitly updated this itera¬ 
tion. Then evaluate / at one point, a feasible stepsize along away from x. 
Using the value of Gjj from the last iteration, and the new function value, com¬ 
pute a second order approximation to gj. Add g 2 to tj 2 . 

Update x B if appropriate. 

Step 5: (Compute mixed second partial derivatives.) Consider each row q, as in Step 3. 
Then for every other q k corresponding to the /th row of G, / < j, evaluate / at a 
feasible point appropriately chosen along q, + q k , and compute the mixed second 
partial derivative G# (this fills in the "pipe" as in Section 1). 

Update x B if appropriate. 

Incorporate all the new second derivatives into the matrix factors of G. 

Note that (n - u)/ir passes through this step, with the same N, will explicitly update all 
of G, using y (n - u) 2 - (n - u) evaluations. In the worst case (u - 0), then, y n(n - 1) 
evaluations are used, on the average. 

The factors of G may be updated at this point using any rank-one or rank-two formula 
which does not force positive definiteness on G on an infinite number of iterations. 

Updating row/column / of G, for / ^ 1, is equivalent to adding a rank-two update of the 
form wej + e t w T , where w /+I * ... - w„ — 0, is the /th column of the identity matrix, so that 
if a scheme utilizing this special structure was devised, at most n updating passes using 
Sorenson’s method [171 would have to be performed. We have not yet devised a procedure 
exploiting this structure. Our code uses a strategy that is reasonably efficient when 
it < — (n — «), and does two rank-one modifications for each row/column updated—the 
second of which wipes out the undesired results of the first. Using the operation count in [181, 
between 2 ir(n - u) 2 + 0(n) and (1 l/3)ir(w - u) 2 + 0(n) operations are required for this 
updating. Full details on updating the Bunch-Parlett factorization are given in [171. 

Step 6: (Check accuracy of derivative approximations.) If G is positive semidefinite and 
the gradient norm rf is sufficiently small, stop; x is an acceptable solution. 
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If G is positive definite but the stepsize is too big relative to the gradient norm. i.e.. 
s > (t>/«K set the failure indicator r = 1 and go to Step 12, where the stepsize will be reduced 
with no movement made. Otherwise, scale the Newton search direction d relative to the ,V~ 
directions. If G is positive definite, go to Step 7, otherwise to Step 8. 

Step 7: (Eigenanalysis when Hessian approximation is not positive definite.) Determine Z 
and A such that Z is orthogonal, A is diagonal, and D = Z\Z T . Let A w denoie 
the most negative eigenvalue of D. If \ qq = 0, D is positive semidefiniie; go to 
Step 8. Otherwise, a direction of strictly negative curvature exists. Solve 
M T y = Z q for it. If the gradient norm is too small, go to Step 9. 

Step 8: (Compute direction of line search.) Solve MiD + ZEZ 0 M T S = —g for 8. where 
£ is a diagonal matrix given by E Jt = (max{|A„l, yin — t/)A,yl - A„] for 

/'= 1,2 . n - w, and A= max IA „ I. £ is thus the "smallest" diagonal 

additions matrix necessary to force positive definiteness on G. If D is positive 
definite, go to Step 10 and search only along the usual Newton direction. If D is 
not positive definite, consider either the "fixed-up Hessian" direction 
1 N*] T d + Qfib or the quadratic curve {.v,l.v, = ,v + tll\*) r d - (?|{,irl + 

t 2 liN + ) T d + where £) n is <he submatrix consisting of the last in — a) rows 

of Q. Note that iN*) T d - Qlg is just the negative gradient in QR coordinates. 
The second vector, ( N*) r d + £?oV. is guaranteed to be a descent direction for / if 
\ qq < 0 and the derivative approximations are sufficiently accurate. Our curve 
search is thus similar to that of McCormick (121 and More and Sorenson [14], 
except that we reverse the t and r 2 terms. Our motivation for doing this is that, 
unlike them, we allow a stepsize t to exceed one. Since the gradient direction is 
most useful only close to the current point, we wanted it to be multiplied by t 
rather than t 2 . 


If the curve is to be searched, go to Step 9; if the "fixed-up" Hessian direction, go to Step 

10. 


Step 9: (Search along a curve.) Use a one variable search method to find a /* ^ r such 
that the point x,. satisfies the two term Armijo-type acceptability condition 
fix,.) - fix) < pit* t> 2 + ( t*) 2 y T g + y(/*) 4 .v T Gy ]. trying t - 1 first. If success¬ 
ful, let x B - x,.. In either case, go to S ep II. 

Step 10:(Regular line search.) Search for a t > 0 such that fix + t{(\'*) r d + @og!) ~ 
fix) < plt(d r Af + 6 r g) + y / 2 8 r C7], trying r- 1 first. If successful, redefine 
x B . 

Step 11:(If the move point is good enough, reduce the stepsize.) If fix BI - 
fix) > -a 2 f) 2 s 2 , there has not been sufficient function value decrease; go to Step 
12. If f(x B ) — fix) < -a 2 fi 2 s 2 , set the local variations failure indicator r - 0. If 
fix b ) - fix) < /3 2 tj 2 , there has been, in addition, sufficient decrease relative to 
the gradient norm, so choose a new smaller stepsize s'and go to Step 13. Other¬ 
wise, x *— x B ; go to Step 1, having moved to x B , but using the same stepsize on 
the next iteration. 
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Step 12:(Failure; try reversing the local variations directions.) Reverse the <r direction 
indicators. If r = 0, this is the first time this step has been encountered with the 
current jcand y, set r = 1 and go to Step 1. Otherwise, if r — 1, let the new step- 
size s' = s/2, set r = 0, x B — x, and go to Step 13. 

Step 13:(Define a new sequence point.) If x ^ x k . k — k + 1, x k — x, s k *— s, x — x B . 
and s — s'. Go to Step 1. 

4. CONVERGENCE 

The asymptotic convergence behavior of our method is similar to that of our modified 
Newton method [10] from which it inherits its coordinate system and searches, and our uncon¬ 
strained method [9] from which it obtains its second derivative matrix approximation approach. 
Proofs for the theorems below are thus obtained by extending results from [9] in the manner of 
[10], so that we omit the complete proofs and indicate only the key points used. 

THEOREM 1: (First order convergence.) If /is continuously differentiable on an open 
convex set containing the bounded set {xl A r x > b , f(x) < /Or 0 )}, then the method con¬ 
verges to a point satisfying the Karush-Kuhn-Tucker first order necessary optimality conditions. 

PROOF. Since fix) is bounded from below, we must have that [s*| is infinite, unless the 
sequence U*} is finite. If so, we can show that the last point found satisfies the optimality 
conditions. Now assume |s*} is infinite. Then (s*) — 0, and, by the acceptability tests at Steps 
6 , 10, and 11, we have that the gradient norm is bounded above by a function of the stepsize, 
so that the projected gradient norm also goes to zero. 

The next two results all require that / be twice continuously differentiable, {**) converge 
to a unique point x, at which strictly complementary slackness holds, and that there is some 
finite bound on the norm of the largest element of the true projected Hessian at 3c, which we 
denote H. 

THEOREM 2: (Second order convergence). If, in addition to the above, V 2 /'(*) satisfies 
a Lipshitz condition in a neighborhood of x , and t is chosen small enough, .v satisfies the 
second order necessary optimality condition that H be positive semidefinite. 

PROOF: This result follows from the explicit updating we tlo of G. If H is not positive 
semidefinite, then U has at least one strictly negative eigenvalue. Since {**) — x. and, by our 
updating, (G*! — H , eventually our curve search direction will be sufficiently close to the 
eigenvector direction associated with that negative eigenvalue, and, if t is small enough, the 
search at Step 9 will be able to be successful infinitely often. Since t is bounded away from 
zero, this would imply that |/(jr*)] — -°o, so that H must be positive semidefinite. The need 
for a "sufficiently small r" is due to the fact that the upper limit of the interval over which t* 
may be chosen to satisfy the acceptability test at Step 9 is a-function of the magnitude of the 
most negative eigenvalue of H , and t must be chosen smaller than that upper limit for the 
proof to follow. 

THEOREM 3: If, in addition to the above, H is sufficiently positive definite in a feasible 
neighborhood of x, the convergence rate of the algorithm is superlinear. If a Lipshitz condition 
also holds, order 2 convergence holds. 
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PROOF: These two results follow from making our approximated Newton line search 
direction computed at Step 8 be sufficiently close to the true direction which would be com¬ 
puted using exact derivatives. The stepsize t = 1 will then be feasible and acceptable, and. by 
McCormick and Ritter [13, Theorem 21, the rates claimed follow. 

5. MODIFYING THE SECOND DERIVATIVE MATRIX 

The algorithm maintains an approximation, G, to the second derivative matrix projected 
onto the coordinates defined by the last in - u ) rows of Q. M and D factors can be main¬ 
tained in the lower triangle of G, and, since G is symmetric, it is stored explicitly in the upper 
triangle of the same array. G can thus be used to aid in adjusting the matrix when the active 
set changes. Two cases must be considered—dropping a constraint from the active set, and 
adding one. 

Dropping a constraint from the active set removes a column from R and a row from A ,+ , 
and adds a row and column to G. An identity column is added to G, requiring no arithmetic. 
(While this may tend to imbalance G, since it will intermix a direction for which no informa¬ 
tion is available with those for which much better approximations are known, we compensate 
by making the new column the first candidate for explicit updating.) The amount of work 
involved in updating Q and R will depend on the number of constraints deleted and their rela¬ 
tive positions in the active set. Removing a single constraint would require from 0 (last con¬ 
straint dropped) to 3un + y u 2 + 0 in) multiplications, 3un + y u 2 + Oin) additions, and 
iu — 1) square roots (first constraint dropped) to update R. and an additional 
3[(« - 2) (w - 1 )/2l + (w - 1)0(1) multiplications/additions for updating Q. using the 
arrangement in [51. 

Adding a constraint to the active set adds a column to R and a row to deflates G from 
u x u to iu - 1) x (h - 1), and modifies Q. n multiplications and n additions are required to 
determine the new column of /?, and then an (n — u) vector must be reduced using orthogonal 
matrices. Using straightforward multiplication, the in - u — 1)2 x 2 Givens' matrices each 
require one square root, 4 multiplications and one addition. Premultiplication of Q involves 
Inin - u) - 3n + 0 in - u ) multiplications and additions. Updating of G is simplified by 
using those transformations that will zero out all but the last entry in the vector, as in Gill and 
Murray [7] We would like to determine triangular factors, LL T , for G, since then premultipli- 
cation by the successive Givens' transformations would change L to lower Hessenberg form, 
and postmultiplication would restore triangularity. 

Now G is recurred as MDM T . with D consisting of 1 x I and 2x2 blocks, and difficulties 
arise when a I x I block is nonposilive or whenever a 2 x 2 block exists, that is, whenever 
there is a nonpositive eigenvalue Fortunately, a fix up of G to be positive definite, so as to 
have LL T factors, is easily done, since whenever G is not positive definite at this step, the com¬ 
plete eigensystem of D has already been determined during the previous execution of Step 7. 
Note that this modification of G would have to be performed once, before each set of con¬ 
straints is added, and it will tend to modify the negative curvature information we try to main¬ 
tain. Again, a least-constrained approach should tend to minimize such modifications. 

The method used follows ideas in Sorenson [171, and is an attempt to replace G with the 
"nearest" positive definite matrix, i.e., G — MDM T is replaced with G — MiD + ZEZ T )M T . 
The exact effect on the second derivative information in G is not clear, but, since we were not 
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able to derive any scheme for preventing the indefiniteness of G which would be more efficient 
than an entire refactorization, we opted for the fix-up approach. Note that this fix up would be 
exactly what we would have obtained if we had used the Bunch-Parlett factorization and forced 
the projected Hessian approximation to be positive definite. 

The fix-up of G and determination of the LL T factor requires between 
y (n — u) (n — u — 1) multiplications, no additions, and (n — u ) square roots in the case of 
all 1 x 1 blocks, to (« - u) 1 + 4(n - u) multiplications, (n - u) 2 + (n - u - 6) additions, 

and (« - u) square roots in the case of all 2 x 2 blocks. The subsequent straightforward mul¬ 
tiplication of the lower triangular factor and its postmultiplication to regain triangularity 

involves 4 (n - u) 2 — (n - u + 2) multiplications and 2(n - u) 2 - (n - u) additions. It 

should be noted that the reduction of the (n - u ) vector 


V| 


to 11v11 u by successive premultiplication by (n — u) 2 x 2 Givens’ matrices of the form 



could be replaced with a single multiplication by the matrix 


C\ 

S| 

0 

0 

0 ... 

S|C 2 

-C|C 2 


0 

0 ... 

S|S 2 C 3 

-ClS 2 C 3 

-c 2 c 3 

■S3 

0 ... 

S|S 2 S 3 C4 

-c,s 2 s 3 c 4 

-c 2 s 3 c 4 

-c 3 c 4 

s 4 ... 


which is a special lower Hessenberg matrix in the sense of 15]. Thus, savings of approximately 
25% of the work involved in premultiplying the lower triangular factor of G could be saved by 
using Lemma IV of Section 2 of 15]. 


If constraints have been added to the active set, N + increases. Due to the complexity of 
updating jV + , especially since several constraints are often added at the same time, N + is 
recomputed at such times. Recomputation of N + requires y u 2 n multiplications and 
y u(u - 1)az additions. 
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6. COMPUTATIONAL RESULTS 

The essential motivating factor in this research was our feeling that a partial explicit 
update of the second derivative matrix might prove to be superior to the usual quasi-Newton 
strategy. We coded our method to see whether numerical experimentation would be consistent 
with that conjecture. (A copy of that code may be obtained from the author.) 

There are actually two different major hypotheses to test here. First, we wanted to see if 
varying w, the number of rows/columns in the update "pipe," could improve convergence 
behavior. We tested this using different values of ir in our new code, and by comparison with 
our older modified Newton method, thus holding the algorithmic structure constant. Second, 
we needed to see whether our new method is competitive with other current software. For this 
comparison, we used a state-of-the-art quasi-Newton method, Murtagh and Saunder's MINOS 
115] algorithm. On a secondary level, we were curious as to how the inclusion of a low rank 
update might affect the results, so we also ran the test problems using Broyden's rank-1 update 
and the complementary DFP rank-2 update. Note that the latter forces positive definiteness on 
G. 


As with similar numerical experimentation with optimization codes, our results cannot be 
considered an absolute conclusion for either hypothesis. The test results using our methods are 
affected by the parameter settings, particularly the starting stepsize. While we tried to find good 
values for each instance, other parameter settings might have been optimal. The path our algo¬ 
rithm took to the solution may differ with it, so that the results may be an artifact of the start¬ 
ing point. Finally, we used the default parameters in MINOS for runs through it. A more 
experienced MINOS user might have been able to speed its convergence by setting its various 
parameters differently. Note that our stopping criterion, since it is a second order one, is more 
stringent than that of MINOS. 

These caveats aside, our experiments do tend to support our conjectures. Table 1 sum¬ 
marizes the number of function evaluations used until convergence (i.e., agreement with the 
optimal / value to 10 decimal places) for the best configuration of our new method, our 
modified Newton method, and MINOS for six popular test problems. The results are 


TABLE 1 — Comparative Function Evaluation Counts 


Problem and 

Jest Form of the New Methoc 

Our Modified 
Newton 
Method 

MINOS 

Source 

IT 

Update 

Evaluations 

Gauthier [3] 
(Colville #7) 

J 

Rank 2 

Ill 

187 

496 

Hydrazine Equil¬ 
ibrium [4] 

7 

None 

213 

223 

561 

Rosenbrock [16] 

2 

None 

109 [9] 

97 

192 

Shell Primal [3] 
(Colville #1) 

0 

None 

44 

50 

60 

Water Quality [10] 

2 

None 

271 

480 

2,292 

Wood [3] 

(Colville #4) 

3 

None 

245 [9] 

447 

585 
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particularly striking for the Water Quality problem, which is a nonconvex 11 variable model 
derived from a model of the Willmette River. The new method performed about the same as 
our modified Newton one on Hydrazine Equilibrium, Rosenbrock's banana function, and Shell 
Primal, and did much better on the other three. 

To further see the effect of modifying it, Figure 1 shows, graphically, the number of 
function evaluations required for convergence on the four linearly constrained problems. Shell 
Primal is a relatively easy cubic problem, and Figure 1A shows that not much benefit results 
from either increasing it or adding a low rank update. The best performance occurred for 
it = 0 and no update; since we start with an identity matrix for G this is nothing more than 
steepest descent. 

Gauthier’s problem is a 16 variable quartic. In Figure IB the convex pattern with respect 
to it that was so marked for unconstrained problems [91 appears for the "no update” and "rank- 
2 update" curves. That is, an intermediate value of it is superior to either the modified Newton 
(it = n) or quasi-Newton (it « 0) strategy. The pattern for the rank-1 update is different, con¬ 
founded in part by local variations moves. Since the method moves to the better of the points 
found by line (or curve) search or points used for derivative approximations, fortuitous values 
of the stepsize and coordinate system sometimes lead to fast convergence. 

Hydrazine equilibrium is a tightly constrained 10 variable problem. The feasible region is 
small, and the method spends much time changing constraint sets (the optimum is strictly inte¬ 
rior to the inequality constraints). Again, particularly lucky searches on local variation moves 
had an impact on the convergence. The pattern of Figure 1C is a bit confusing, although the 
modified Newton-like strategy of larger it seems to be desirable. 

Finally, the performance of the Water Quality problem is illustrated in Figure ID. Since a 
rank-1 or rank-2 update cannot be performed unless the same constraint set is used for two 
consecutive iterations, the constant changes of constraint sets meant that such updates were 
never performed. Here the "U" shaped patterns appear again, lending support to the idea that 
some explicit updating is better than none, but that a complete update is usually not 
worthwhile. This is also a rather hard problem; note that MINOS required 2,292 function 
evaluations to converge. 

Table 2 shows the relative percentage of time used by each part of our algorithm. Our 
code is an experimental one, written by the author, so that a more skillful implementation 
might redistribute these somewhat. Test runs were performed on the University of Pittsburgh’s 
DEC 1099 in the usual batch environment. In this situation, times tend to vary by ±15%. 
The standardized times reported in the last two columns were computed relative to the 28.48 
seconds required for Colville’s [3] timing routine (average of 3 runs). 

The bulk of the time required by the new method was for the numerical derivative 
approximations. This was reassuring, since a basic idea in using a nonderivative method, such 
as ours, is to trade off increased computer time for the human time necessary to compute the deriva¬ 
tives analytically. 

Finally, the standardized time ratio for our new method compares favorably with that of 
MINOS. The percentage reduction in time for our method, relative to MINOS, ranges from 
37-70% for the constrained problems and is about 90% for the unconstrained ones. Considered 
together with the additional time necessary to derive and program the derivatives for a code 
such as MINOS, we can conservatively conclude that our new method is at least competitive 
with a representative, state-of-the-art, quasi-Newton algorithm. 
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TABLE 2 — Timing Results 

Time <sec ) 


Setup^Findingj 
Feasible Poinl 


|Consiraint Handling 


Standardized 

Time 

Ratio (total) 


Standardized 

Time 

Ratio (total) 


Gauthier 
Hydrazine 
Rosen brock 
Shell Primal 
Water Quality) 
Wood 


00557 

0.0331 

0.0015(19! 

0.0109 

0.1131 

0.0066(19] 


0.0885 

0.0572 

0.0330 

0.0372 

0.1921 


7. CONCLUSIONS 

Our thesis, in this paper, is that explicitly updating a few columns of the approximated 
(projected) matrix of second partial derivatives at each iteration may yield convergence 
behavior superior to that of traditional quasi-Newton algorithms. We presented a method for 
linearly constrained optimization, incorporating that idea into our existing algorithmic frame¬ 
work. The user specifies the number of rows/column to be updated at each iteration. Theoret¬ 
ical convergence proofs show that such a method retains the properties of our modified Newton 
algorithm—convergence to a second-order Kuhn-Tucker point, and superlinear rate. Numerical 
experiments indicate that this approach is competitive with a state-of-the-art quasi-Newton 
code, and that a small number of explicitly updated rows/columns (usually 2 or 3) can 
significantly improve the speed of convergence. 
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ABSTRACT 

In this paper, we develop efficient interactive methods for the solution of 
bicriteria nonlinear programming problems. The methods do not require 
trade-off information from the decision maker, pose less cognitive burden and 
converge to the "best compromise solution" fast. Two methods, called the 
paired comparison method and comparative trade-ofT method, are presented 
with examples. A real application of the interactive method to a bicriteria prob¬ 
lem that arose in the planning of the cardiovascular disease control program in 
the U.S. Air Force is also presented. 


INTRODUCTION 

Bicriteria mathematical programming (BCMP) problems as a first generalization of single 
criterion nonlinear programs have been studied by many researchers [1, 2, 3, 4, 7, 13, 17J. 
The small dimensionality of the problem permits a visual presentation of the "payoff set," which 
partly explains such widespread interest. Further, when we restrict ourselves to "efficient solu¬ 
tions," the complementary roles of the two criteria yields additional computational power. 

The BCMP problem may be stated as follows: 

(1) VMAX {/,(*)./ 2 0c» 

Subject to g,(x) <0 i - 1. m 

where x is an n dimensional vector of decision variables, ,/j and f 2 are real valued criterion func¬ 
tions and g ,'s represent a set of nonlinear constraints. Let 5 - (x|&(x) < 0) denote the feasi¬ 
ble region. An optimal solution of BCMP will be taken to be a "best compromise solution" that 
maximizes the preferences of a Decision Maker (DM); that is, a feasible solution that solves 
the following program. 
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(2) Max u\f\(x).f 2 (x)) 

Subject to x € S 

where U is the DM's preference function or value function defined over the criterion values 
</j, f 2 ) such that for any two alternatives *i and x 2 > [/i (jc, ), / 2 (at, )] > U[/|(x 2 ), / 2 (jc 2 )I. if 

the DM prefers x t to x 2 . 

In one of the earlier papers addressed to the bicriteria problem. Geoffrion [7] considered 
the "scalar maximization approach,” i.e., maximizing a convex combination of the two criteria. 
Me showed how such problems can be numerically solved using parametric programming algo¬ 
rithms. Bacopoulas and Singer (21 on the other hand, used the "constrained criteria approach." 
i.e., maximizing one criterion keeping the other criterion as a constraint. They showed how all 
the efficient solutions can be generated by parametrically varying the level of the constrained 
criterion over a particular interval. Some further refinements and alternate proofs appear in 
Gearhart [6] and Benson [31. Pasternak and Passy [13], addressed the special case when the 
variables are zero-one and devised a special enumerative algorithm for generating the efficient 
solutions. Choo and Atkins [4] address the case when the criteria are linear fractional func¬ 
tions. Adulbhan [1] considers the special case when the criteria are linear. In all these 
methods, the emphasis is on generating all the efficient solutions of BCMP. But Walker 117] 
considers the "interactive" approach and uses the Generalized Lagrange Multiplier method of 
Everett [5], Walker’s method is slow in convergence and seeks difficult interaction from the 
DM. 


In this paper we will consider interactive approaches to the solution of BCMP where the 
preference function U is only implicitly known. Such an approach in the more general context 
of multiple criteria problems have been addressed by Geoffrion, Dyer and Feinberg [8] and 
Zionts [19], However, Wallenius [18] has observed that such methods pose considerable cogni¬ 
tive burden on the part of the DM; also the procedures converge to the optimal solution rather 
slowly. In this paper, we develop interactive procedures that pose less cognitive burden to the 
DM and converge to the optimal solution fast. 

Section 2 presents the principal results that form the basis of the interactive methods 
developed for BCMP. Section 3 discusses one of the interactive approaches called the Paired 
Comparison Method. Section 4 discusses another approach called the Comparative Trade off 
Method. Section 5 discusses a real application of the interactive method to a bicriteria problem 
that arose in the control of cardiovascular disease in the U.S. Air Force. 

2. PRINCIPAL RESULTS 

In this section, we will establish a number of results that form the basis of the new 
interactive procedures which are develolped in the later sections. We will assume throughout 
that the feasible region S is a compact convex set; the criterion functions f\ and f 2 are 
differentiable concave functions and the preference function U is differentiable, increasing and 
strongly quasiconcave defined in the following manner. 

DEFINITION 1: f:S —> R\ where S is a nonempty convex set in R" is said to be 
strongly quasiconcave on S, if for each x t . x 2 € 5 with X\ ^ x 2 and for every A € (0, 1), we 
have 

(3) /[AX| + (1 - X)x 2 ] > Min I./‘(jci), fix 2 )\. 
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The convexity assumptions on criteria and feasible region are made only to allow the 
methods of convex programming to be applicable. The quasiconcavity assumption of U is well 
supported by preference theory (see Intriligator [10]). 

Use of the concavity assumptions guarantee the attainment of maxima and quasiconcavity 
ensures us that local optima are also global optima. In addition it is sufficient to restrict 
ourselves to efficient solutions defined below as candidates for the optima. 

DEFINITION 2: A solution x° € 5 is said to be efficient to BCMP. if /|(.v) > /, (V) for 
some x € S => / 2 (x°> > / 2 (x). 

We will now develop procedures for generating the efficient solutions and then a method 
of "efficiently" searching among the efficient solutions for the "optimal" solution. 

(4) Define the payoff set Y - {.v|/(x) — y for x € 5}. 

(5) Let A = |y|/(x) ^ .v for x € 5). Obviously, Y C A. 

LEMMA 1: If Sis a convex set and /is concave, then A will be a convex set (12]. 

Consider the following single objective nonlinear programs: 

PI: Max/, Or) P2: Max / 2 (x) 

Subject to: x 6 S Subject to: x € S 

Let the optimal values of (PI) and (P2) be v*and w* respectively. Now consider the fol¬ 
lowing mathematical programs: 

P v : Max/ 2 (x) (?„.: Max ,/j(x) 

Subject to: x € S Subject to: x € S 

/, (X) > v /,(x) > H' 

Let x* solve Q w for w = w*. Then /, (xj) is the minimum achievable value for ,/j 
without sacrificing any achievement on / 2 , while v*is the maximum value of ,/j at the expense 
of / 2 . Hence, the range of achievable values for /|, denoted by v, is given by (,/j(x 2 ), v*] and 
the best compromise value for /, lies in this range. 

THEOREM 1: In the optimal solution to the program P- where v € l/j(xj), v*), the 
f constraint /,(x) > v will be a binding constraint, i.e., if xsolves P x then /,(jf) - v. 

PROOF. Assume the contrary, i.e.. 

(6) /,(x) > v > /,(x?) 

w* being the absolute maximum of f 2 over S, 

(7) w* - / 2 (xf) > / 2 (x) 
we will consider the two cases 

(8) (i) Let / 2 (xJ) - / 2 (x) 

(6) and (7) taken together contradict the fact that xj solves Qt. 
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(9) (ii) Let / 2 (xJ) > / 2 (x) 

consider the line segment [x, x*l which lies in the convex set S. Over this line segment, the 
concave functions J\ and f 2 are unimodal. Because of "strict” inequality in (6) and (9). there 
exists a feasible solution x' on the line segment (x, xj] such that 

(10) /,(.v) > /,(x') > v 
,/j(x) < /j(x’) < / 2 (xf). 

In equality (10) indicates that v' is a feasible solution to P f with / 2 (x') > / 2 (x) contradicting 
the optimality of x. 

Hence. /|(x) = v. 

THEOREM 2: x 6 5 is efficient if and only if x solves P s where v € [j\ (x‘). v*]. 

PROOF 

Sufficiency. Let x solve where v € lf\ix\). v*J. We will show that x is efficient. Con¬ 
sider an x € S such that 

f\ (x) > /|(x) - v. 

It is clear that x is feasible to / > 5 . But x cannot be optimal to P x : as per theorem I every 
optimal solution x° satifies /|(x°) - v. Hence, f 2 (x) < f 2 (x). Consider an x € 5 such that 
/ 2 (x) > f 2 (x). 

Since x is optimal to P xcannot be feasible to P K . Hence, /|(x) < J\ix). Consequently, x 
must be efficient. 

Necessity: Let £be the set of all efficient solutions. Suppose x € £. Obviously /,(.v) ^ 
./,(xj) as, otherwise, xj will dominate x. Let f\(x) - v. Assume that xdoes not solve P but 
some other x' € 5 solves P v By theorem 1. /](x') - v - /,(x) since x' is optimal for P 
f 2 (x‘) ^ f 2 (x): but x € £ —> /j(x') < f 2 (x) and hence /j(x') - / 2 (x). Consequently, x 
solves £-. 

Theorem 2 enables us to generate the entire set of efficient solutions by parametrically 
solving P v By theorem 1, the generated solutions will have specific levels of attainment of 
namely v. The next theorem would give some unimodality property enabling us to make an 
efficient parametrization. 

Consider the problem where we maximize V over the pay-off set K keeping f, fixed at v; 


P3: Max U(x,f 2 ) 

«v./ 2 >« y 

where (v,/ 2 ) € Tis a subset of Tdefined by Uv,/ 2 )|/ 2 - / 2 (x) for x € Sand f\ix) - v| 
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By Theorem 1 and the fact that U is increasing in f h it is easy to see that an optimal solu¬ 
tion to (P3) is also an optimal solution to program <P K ) where v € (.vj). v*]. We now estab¬ 

lish the unimodality of defined below: 

THEOREM 3: Let V be a convex subset of R 2 and is a strongly quasiconcave 

increasing function of J) and J\ defined over Y'. then 
,?(v) = Max t/fv./j) 

is strongly quasiconcave, i.e., unimodal in v. 

PROOF: Let 

#(v,) = Max U(v } ./ 2 ) - U(\ j,w 3 ) (say) 
j»(v 4 ) = ^ Max ^ l/(v 4 ,/ 2 ) = L'(v 4 ,w 4 ) (say). 

Assuming boundedness of functions the maxima are attained; by convexity of ) and strong 
quasiconcavity of U they are attained uniquely. 

Note that (v 3 ,w 3 ) 6 Y' and (v 4 ,w 4 ) € Y' by convexity of Y\ 

[Av 3 -Ml - A)v 4 , Ate, MI - A)h' 4 ) also € Y. for 0 < A ^ 1. 
j?[v 3 + (1 - A)v 4 ] - Max U [A v 3 + (1 - A)v 4> / 2 1 

(By convexity of Y) ^ t/(Av, + (1 - A)v 4 .Aw 3 + (1 - A)h- 4 1 

(By strong quasi- > Min [LMvj.w,,). f/(v 4 .w 4 )l 
concavity of U) = Min lj?(v 3 ), ^(v 4 )] 

Hence, j?(v) is unimodal in v. 

In the proof of the above theorem, we have heavily exploited the convexity of > When 
we want to apply the result to K the pay-off set, we find that the assumption of convexity of 
pay-off set may be overly restrictive. For example our assumptions of convexity of S and con¬ 
cavity of //s are not sufficient to guarantee convexity of Y. as the following simple example 
would illustrate. 

EXAMPLE: S - x > 0: /- (/,,/j) - (x -log x) Obviously, Sis a convex set and s 
are concave. But the pay-off set Y is the graph of -log x which obviously is not a convex set 
[ 121 . 

However, we will presently establish the fact that we can extend the pay-off set Y to the set .4 
defined in Equation (5) without affecting the results. Lemma 1 established the convexity of 
the set A. Hence identifying A with Y\ the result of theorem 3 holds even when the pay-off 
set is nonconvex. 
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We will establish in the next theorem, the equivalence of solutions to the following two 
problems. 

P4: Max y U(a./ 2 ) P5: ( Max^ U(a,f 2 ) 


THEOREM 4: The optimal solutions to (PA) and (PS) are identical, i.e., if Y m be the set 
of optimal solutions to (P4) and A m be the set of optimal solutions to (P5). then Y m = A m . 

PROOF: 

(a) We will first prove that every element of A m is an element of Y m i.e., A„ D Y m . 
Let y* = ( y*,y 2 ) maximize P5 for a - yf, i.e., y* € A m . We have to show that 
.v* € Y and hence y* € Y m . Since A m C A , by definition 3 and x*, such that 
f(x*) - (ft.JV ^ y* Note that y* € Y is fix*) — y*. Assume that 
fix*) > y*, i.e.. 

Ml) ft>yt.ft>yt( say). 

Since U(f\.f 2 ) is increasing in f\ and / 2 , (11) gives 

(12) (/(yf-yj) < U(ylff). 

But from the fact that y* maximizes U (y*,/ 2 ) over € A we get 

(13) Uiylyp > Uiylfi). 

(13) contradicts (12); hence our assumption (11) must be invalid and fix*) - y*. By 
definition of the set K. y* 6 Y. Being a maximizer of (P5) for a = yf. y* € Y m and hence 
A m D Y n . 

(b) Assume y* € Y m . We will show that y* € A m as well or Y m C A„. Since 
Y„, C K c A, we get y* € A. Assume y* is not optimal to (P5). Hence 3 
another y' = (y*,y 2 ) such that 

Viy'.yP - ( Max^ U ly*,/ 2 ] 

(14) U (y*.y:’) > U(y’,f 2 ) for any (yf./ 2 ) € A 

and in particular > U (y’.yj). 

Now (yf,yj) € A m : by part (a) of the proof we will have (y*. Y' 2 ) € Y. By optimality of 
(yf.yp for (P4) 

(15) U (y*,y 2 ') < U (y*.y 2 ) 

(14) and (15) imply that 

Uiy*.y 2 ) = V(y*,y 2 ) 

or in other words (yf.y?) optimizes (P5) as well. Hence 
y* 6 A m or Y m C A m . 
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THEOREM 5: With the assumptions mentioned earlier the following relationship is true: 

Max Li (v,./' 2 ) = .g(v) = Max U(v,f 2 ) 

Subject to: ,/j(.v) > v 
.v € S 

PROOF: Follows directly from Theorems 3 and 4. 

3. PAIRED COMPARISON METHOD 

Using the results of Section 2 we will develop a procedure needing only a /wired compari¬ 
son from the DM. i.e., given any two feasible solutions and their outcomes. sa\ v' 11 and >•*-' 
the DM only has to specify whether 
(16) /" > y or y y >■'" or both 

where denotes "preferred to.” 

Interaction of this form is presumably much less demanding than asking the DM to 
specify his local tradeoffs. We will presently indicate how (16) can be used to eliminate a por¬ 
tion of the efficient set and progressively converge to the "best compromise solution." 

Based on the theorems given in Section 2, the BCMP problem has been reduced to deter¬ 
mining the maximum of #(v) where v belongs to the interval [v,. v*l where v, = ./j(.v*>. How¬ 
ever, #(v) is not known explicitly since U is not known. But, using a search technique which 
requires only functional comparison and not function values, we can still solve the BCMP prob¬ 
lem, using the following region elimination concept: 

For a unimodal function j?(v) defined over a finite interval (v,.v*), let v 4 and v fl be two 
points in the interval such that v 4 < v B . Then, g(v 4 ) < g(v 8 ) implies that the maximum of 
g(\) will not lie in the interval (v,,v 4 ). On the other hand, ,g(v 4 ) > g(v B ) implies that the 
maximum will not lie in the interval (v B ,v*). 

General Steps of the Algorithm 

Step 1: Solve PI: Max ./j(.v), subject to x 6 S. 

Set max f\(x) = v* 

Solve P2: Max ./' 2 (jc), subject to x € S. 

Set max / 2 (x) = w* 

Step 2: Solve Q*:. Max ./j(x), subject to x € 5and / 2 (x) > w*. 

Set the maximum of f\ v ; . Now the optimal value 
of ,/j to the BCMP lies between v, and v*. 

Step 3: Choose two values, v,, and v s , such that v, < x A < v B < v*. Solve the problem 
P y : Maximize / 2 (x), subject to x € S and ,/j(x) ^ v for v = x 4 and v = v fl . let 
glv) - Max / 2 (x) for Py. 
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Step 4: Let y u = [v^.gfvy], y <2) = [v B ,g(v B )]. Given V n and y*’, the DM is asked to 
specify whether y iU >. (preferred to) v* 21 , or j' ,2) ^ > ,ni , or indifferent. Usine the 
DM’s response, a portion of the efficient set can be eliminated. 

If y n > y (2 \ then set v„ =» v B ; go to step 5. 

If y 21 V y n , then set v, = v 4 ; go to step 5. 

Step 5: If Ivj - vj < e (chosen small number) stop; otherwise return to step 3. 

Figure 1 illustrates the paired comparison method. Here, y <21 > .v ( " since U(y t2 ') > 
U(y {U ). Hence, the maximum of U cannot lie in the interval [v/. v j. Eliminating the region 
(v/.vy, the interval of uncertainty where the optimum lies reduces to (v^.v*). If Golden Sec¬ 
tion Section Search is used to generate the two new points between v 4 and v* one of the new 
points will turn out to be v B . With Golden Section Search, we will be able to bracket the best 
compromise solution y° to less than 10% of the original interval on v with just 5 paired compar¬ 
isons from the DM; in 10 paired comparisons, the optimal solution will be within 1% of the ori¬ 
ginal interval of uncertainty [v ( , v*]. Thus, for a specified level of uncertainty we have a finite, 
rapidly convergent procedure using only paired comparison of two efficient solutions. 
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Numerical Example 

Consider the following bicriteria program [7] with 

/,0c) = 32 -40 x 2 + 23 x 3 — 7 x 4 

/ 2 (.x) = 32 — 10 X| — 4 x 2 + 7 X} — 7 x 4 
and the feasible region S is specified by the following linear equality constraints in detached 
coefficient form in addition to the usual nonnegativity restrictions: 



*2 

x} 

*4 

*5 

*6 

b 

I 

0 

-23/50 

27/50 

0 

0 

132/50 

0 

1 

-65/100 

35/100 

0 

0 

60/100 

0 

0 

1/25 

1/25 

1 

0 

16/25 

0 

0 

15/50 

-35/50 

0 

1 

40/50 


In order to simulate the interaction process, we will assume that the DM's implicit prefer¬ 
ence function (which guides him in the selection process) is the following quasiconcave func¬ 
tion 

V(fuh) - fl n h- 

The true optimal solution is x* = 

= (1.2857. 0. 0.7895. 3.1805, 0.4812. 2.7895) with /* = (ft, ft) 

= (27.895, 2.406) and U(ft, ft) - 22.13. 

Solution Using Paired Comparison Method 

Step 1: Solving (PI) and (P2) we get v* = 60 and w* = 3.2 
Hence v u = 60 

Step 2: Solving (£);.) yields v, = 8 
Iteration 1 

Step 3: Using Golden Section ratios 0.618 and 0.382, v A = 27.864 and v 8 = 40.136. Solve 
the problem P y for v = v< and v fl . g(v 4 ) = 2.407 and g(v fi ) = 1.707 

Step 4: Interact with the DM and seek his paired comparison between y il) = (27.864, 
2.407) and y l2) = (40.136, 1.707). Assuming that DM is guided by the implicit 
preference function U(f.fi) = f? 3 / 2 , £/[y <,, l = 22.12 and f/(y (2) ] = 20.01. By 
the property of the preference function, /" > y i2) if t/iy 1 ’] > £/ly <2 ’]. Hence, 
would be preferred to y <2> Consequently, v B will be updated to v s - 40.136, 
thereby eliminating a portion of the efficient frontier. 

Step 5: Assuming that the interval of uncertainty is to be reduced to 10%, the search has 
to be continued further on the reduced efficient frontier. For brevity of presenta¬ 
tion the results of further calculations are summarized in Table 1. 
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TABLE 1 — Solution of Bicriteria Problem Using 
Paired Comparison Method 


Iteration # 




V B 



1/ [>■'"] 

t/|.v'-”l 

Remaining 
Interval 
('•/„ of Orig- 

0 

8 

60 

_ 

_ 

_ 

_ 

_ 

- 

100 

' 

8 

40.136 

27.864 

40.136 

(27.864, 

2.407) 

(40.136. 

1.707) 

22.12 

20.01 

61.8 

2 

20.276 

40.136 

20.276 

27 864 

(20.276. 

2.841) 

(27.864, 

2.407) 

21.12 

22.12 

38.2 

3 

20.276 

32.549 

27.864 

32.549 

(27.864. 

2.407) 

(32.549. 

2.140) 

22.12 

21.82 

23.6 

4 

24.964 

32.549 

24.964 

27.864 

(24.964. 

2.573) 

<27.864, 

2.407) 

21.97 

22.12 

14.6 

5 

24.964 

29.652 

27.864 

29.652 

(27.864. 

2.407) 

(29.652. 

| 2.306) 

22.12 

22.09 

9.0 


4. COMPARATIVE TRADEOFF METHOD 

The local trade-off at a feasible point x € S, between criteria 2 and 1 is defined to be 



evaluated at x. 

T 2 1 measures the loss in criterion f 2 from the current level the DM is prepared to trade¬ 
off for unit gain in criterion f x . 

If the DM can provide local trade-offs then the problem can be solved interactively using 
Geoffrion-Dyer-Feinberg [8] type approaches. However, it has been repeatedly observed in 
practice that precise tradeoffs are difficult to provide. However, the DM can provide "impre¬ 
cise" tradeoffs much more easily, e.g., "interval estimates" in place of "point estimates." For the 
bicriteria case, an estimate of the local tradeoff that is easier than the interval estimate is the 
"comparative estimate." The DM is provided with a number and will be asked to respond 
whether or not he would be prepared to tradeoff more than, less than, or equal to the specified 
number in criterion J’ 2 for a unit gain in criterion f\. Based on the DM’s response, a portion of 
the efficient frontier can be eliminated. We will presently indicate how the results of Section 2 
can be used to provide the basis for the "comparative trade-off method." 

At any feasible point x € S, let X 2 t be the perturbation in f 2 from its current level per 
unit perturbation in /,. Limiting ourselves only to efficient solutions, we find that Theorem 2 
provides us a method of generating efficient solutions; also Theorem 1 specifies that the con¬ 
straint /|(x) > v in program P- will be binding in every optimal solution. Hence, x 2) is 
obtainable as the negative of the Lagrange multiplier associated with the constraint f x (x) > v. 
Note that X 2I is part of the solution output of most mathematical programming algorithms. 
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X 2t has the following interpretation: X 2 | measures the loss in / 2 from its current level that 
must be suffered per unit gain in /,. On the other hand T 2X measures the loss in f 2 the DM is 
prepared to trade-off per unit gain in f\. Naturally, these two quantities must be related to the 
optimality of BCMP and the precise relation is the subject of our next theorem. 

THEOREM 6: 

Let x* be the optimal solution to BCMP with fix*) = (/*. /*). Let x be any efficient 
solution. Then 

(a) x* is optimal to BCMP if and only if \ 21 = T 2 |, at x*. 

(b) At x if k 2 \ < T 2 i then f\ < ./?; similarly, if X 2 | > T 2 then /, > f\. 

PROOF: Part (a) Sufficiency: At x* let X 2 | = T 2i . By definition of r 21 , at x* the DM is 
indifferent to unit gain in ,/j and A f 2 units of loss in / 2 where A f 2 = J'f^~Jy/ ln 

other words 

(18) t/i(x*>,/ 2 (x*)| ~ [/,(x‘) + 8|,/ 2 (x*) - A J 2 ■ 8,) 

where denotes the indifference of the DM and 8i is the differential change in ,/j. When 
A 2 i = Tji, in the neighborhood of x*. a differential change 8, in ,/j is accompanied by a 
differential change of -(A f 2 ) • 8| in / 2 , as per the definition of X 2I . Hence, points in the 
neighborhood of x* are equally preferred to x*. hence, x* is the local optimum. By quasicon¬ 
cavity of U, x* must be the global optimum too. 

Necessity: Let x* be optimal to BCMP but X 2 i ^ r 2) say X 2I < 7" 2 |. Around a neighbor¬ 
hood of x*, 3 a point l./j(x*) + 8,, / 2 (x*) - 8|X 2 |1 where 8| is a differential change in ,/j as 
per the definition of X 2i . At X 2 | < 7^, 

(19) (/ 2 (x*) - 8|X 21 ) > (/ 2 (x*) - 8, r 2l ) for 8, > 0. 

Since V is strictly increasing in both f\ and / 2 , inequalities (18) and (19) taken together imply 
that 

(20) 1/ 1 (x*) + 8|,/ 2 (x*) - 8|X 21 ] > l/,(x*) + 8|,/ 2 (x*) - S,r 2 |] 

— I/,(x‘)./ 2 (x*)J. 

By transitivity of preferences 

(21) L/j(x‘) + 8„ / 2 (x‘) - 8|X 21 ) > l/,(x*)./ 2 (x*)l. 

Obviously, (21) contradicts the optimality of x*. Hence, X 2 | < T 2l . Similarly we can show that 
X 2 | > r 2 |. Hence the proof. 

Part (b) Let X 2 | < r 2) . By (21) an increase in /, from its current value leads to a pre¬ 
ferred solution. Hence, an increase in /j is a locally increasing direction for U along the 
efficient frontier. Also x being an efficient solution is a solution to program (P v ) by Theorem 2 
and f\ (x) - v by Theorem 1. We noted earlier that x is a solution to program (P3) also, i.e.. 
Max t/(v,/ 2 ). Theorem 5 established the unimodality of g(v) which is the "value" of 
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optimal solution to program (P3); as f x (x ) - v, the value of optimal solution to program (P3) 
is unimodal in f x (x). By the property of unimodality, U increases along a locally improving 
direction continuously up to the optimal solution and continuously decreases in the other direc¬ 
tion. Hence, /j(x*) > /j(x) or/f > /,. Similarly, we can prove the other case also. 

Note that Theorem 6 provides a convenient way of ruling out a portion of the efficient 
frontier. All that is needed is an interaction with the DM whether or not X 2 | < T 2 \- Depending 
on his response we can narrow down the search of efficient points to only those with /, value 
greater than or less than the current /, value. We wish to point out here that the Surrogate 
Worth Tradeoff method due to Haimes, Hall and Freedman [9) is essentially based on a result 
similar to part (a) of Theorem 6. However, they do not have a result comparable to part (b) of 
the Theorem. Thus no "region elimination" is possible in their method. Also the DM has to 
specify how far T 2X is greater than or less than A 2 | on a subjective scale of + 10 to -10. 


We now formally state the procedure as follows: 

Step 1: Solve (PI) and (P2) and determine v* and w*. 

Step 2: Solve (Q w .) and determine / t (xJ). Set v, - /j(xj). 

Step 3: Solve P v for v - v^ where v t < v A < v„. Determine X 2 i, ie., the Lagrange mul¬ 
tiplier corresponding to the constraint /j(x) > v A in the program P v . 

Step 4: Interact with the DM and present him with A 21 and seek his comparative trade-off 
7*21* i.e., whether T 2 \ £ X 2 |. Depending on the outcome of interaction eliminate a 
portion of the efficient frontier using the following rule: 

If X 2( < f 2 „ update v, - v,,; go to step 5. 

If X 2 | > f 2 |, update v„ - v^; go to step 5. 

If X 2 j - r 2 |, stop; the current solution to P v is optimal to BCMP. 

Step 5: If |v, - vj < ( (chosen small number) stop; else go to step 3. 

Note that selection of v^ values can be done efficiently using the midpoint of the interval. 
For example, with S stages of interaction, we will be able to bracket the optimal solution to less 
than 4% of the original interval. In 7 stages, we will be within 1%, etc. Thus we have a finite, 
rapidly converging procedure using only "comparative tradeoffs." We shall illustrate the method 
using the sdme example given in Section 3. 

Numerical Example 

Solution Using Comparative Trade-Off Method: 

Step 1: As before v* - 60; w* - 3.2. Hence v„ •« 60. 

Step 2: As before v, - 8. 
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Iteration 1 

Step 3: Using the 'Bisection’ method (i.e., midpoint values) to choose v^ values, v A = 34. 
x A which solves P v for v - v^ has / values of (34,2.057); also the Lagrange multi¬ 
plier A. 21 = 0.0571. 

Step 4. Interact with the DM and seek his comparative estimate of local trade-off T 2 \, i.e., 
whether T 2i > A 2) , T 2I < A 2 , or T 2 , - A 2I . Assuming that the DM is guided by 
the implicit preference function 

U(f } .f 2 ) - /, 2/3 / 2 ; we get 

~j- = 2/3 /f' 1/3 f 2 and — = /, 2/3 . Evaluating at x A , 

~ - j (34)- ,/3 (2.057) - 0.4233 and 

f£-34*- .0.4*1. Hence. -^--0.0403. 

Thus, A 2 1 «■ 0.0571 > r 2) - 0.0403. Hence, update v„ to \ A - 34 

Step 5: As the limits v, and v„ are not close enough, the search has to be continued 
further on the reduced efficient frontier. The details of the calculation for further 
stages are summarized in Table 2. 


TABLE 2 — Solution of Bicriteria Problem Using 
Comparative Tradeoff Method 


Iteration # f u, [ 

v. 

1 

*2! 

r 2 . 

Remain Interval J 
(% of original interval) 1 

1 0 1 8 ! 

60 


_ 

- 

100 I 

! 1 I 8 ! 

34 

1 34 (34,2.057) 

0.0571 

00403 

50 

• 2 21 

34 

21 (21.2.800) 

0.0571 

0.0889 

25 

1 3 | 27.5 | 

34 

| 27.5 j (27.5.2.429) 

00571 

0.0589 

12.5 

| 4 j 27.5 | 

30.75 

| 30.75 | (30 74,2.243) | 

0.0571 

0.0486 

6.25 | 


5. BICRITERIA OPTIMIZATION APPLIED TO THE U.S. AIR FORCE 
CARDIOVASCULAR DISEASE CONTROL PROGRAM 

In this section we will demonstrate an application of the paired comparison method to a 
bicriteria optimization problem of designing a cardiovascular disease control program for the 
U.S. Air Force personnel, this part of a study currently undertaken by Purdue’s School of 
Industrial Engineering, under contract with the USAF School of Aerospace Medicine. 

The United States Air Force (USAF) is planning a comprehensive program for reducing 
the incidence of cardiovascular disease (CVD) among its active duty personnel. This program, 
known as the Health Evaluation and Risk Tabulation (HEART) program, would consist of 
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(i) an extensive educational program about CVD, its causes and prevention, 

(ii) a risk screening program to identify those AF personnel who have a high susceptibil¬ 
ity for developing CVD, and 

(iii) a risk intervention program to treat such high risk personnel and reduce their CVD 
risk by medications and/or behavior modifications. 

In order for such a program to be viable, the expected benefits resulting from reduced 
incidence of heart disease must outweigh the cost of the HEART program. 

Risk Identification: A general cardiovascular risk profile, based on the Framingham Study 
[11], is presently being considered as the means of identifying high-risk personnel for treat¬ 
ment. With this approach a logistic regression model would be applied to each individual and 
would yield his CVD risk; namely, the probability of manifestation of cardiovascular disease in 
the individual in an eight-year time period. 


Risk Intervention: By rank ordering the CVD risk scores of USAF personnel, a certain 
percentage of the highest risk individuals are selected for specialized face to face therapy for 
modifying their elevated values of risk factors. At present, those high-risk individuals with 
elevated systolic blood pressures will be treated with medication, while those with elevated 
cholesterol levels will be treated with diet modification only. Both individual and group thera¬ 
pies will be used to persuade smokers to quit smoking. 

The recurring cost of the HEART program, excluding the initial start-up costs, is 
estimated between 7 and 9 million dollars a year depending on the total size of the therapy 
group. Hence, before implementing the HEART program on all airforce bases, the Air Force 
was interested in answers to the following question: 


Given the total size of the therapy group, determine the optimal 
number of Air Force personnel to be selected in each age group 
of flyers and nonflyers such that the total "cost-effectiveness" of 
the HEART program is maximized. 


The therapy selection problem has indeed two competing objectives. One objective is to 
minimize the USAF CVD risk and the other objective is to minimize the HEART program 
cost. By increasing the size of the therapy program, the CVD incidence in the Air Force will be 
progressively reduced but this can be achieved only at the cost of increased HEART program 
expenditure. If we follow a policy of selecting the high risk personnel from each age group for 
risk intervention, the graph of CVD incidence versus therapy size will be a monotonically 
decreasing function. Hence, there is an optimal expenditure level for the HEART program, an 
optimal level of CVD incidence, an optimal therapy size and an optimal selection strategy for 
that therapy size. Determination of all these quantities is a nontrivial problem for which we 
developed the following bicriteria mathematical programming model which was solved using the 
paired comparison method developed in Section 3. 
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BCMP Model 

Minimize 
Minimize 
Subject to 


Where: 

x — selection strategy (x\,xi .x 17 ) for each age group of flyers and nonflyers. 

f\(x) —USAF CVD risk for a given selection strategy x 

fi(x) —HEART program cost for a given selection strategy x 

S = Maximum size of the therapy group (fixed) 

a, - Minimum number to be selected from age group j (flyers and nonflyers); 

7-1.8 correspond to flyers in age groups 20-24.55-59, respecitvely; 

7-11 .17 correspond to nonflyers in age groups 15-19.55-59, respectively. 

bj - Maximum number that can be selected from age group j. 

Even though f}(x) could be determined analytically, / t (x) is indeed a random variable 
because of the unpredictability of the effectiveness of the risk intervention programs on modi¬ 
fying the risk factor levels. Hence, a comprehensive simulation model known as P-HEART 
(Purdue Health Evaluation and Risk Tabulation) was developed by the authors [14]. P- 
HEART simulates the USAF population and performs risk identification, therapy selection and 
risk intervention. The description of the simulation model, its validation and findings are given 
in Reference [14], Through extensive simulation runs, the expected value of f\(x) for different 
selection strategies were estimated and a polynomial equation was fitted to get analytical forms 
for f x (x) by age group and flying status. 

The paired comparison method was then programmed within the frame work of interac¬ 
tive solution of BCMP for minimizing CVD risk and HEART program cost. The Generalized 
Reduced Gradient Algorithm was used to solve the nonlinear programs in Steps 1, 2, and 3 of 
the interactive algorithm. The interactive program was extensively tested at Purdue using Air 
Force Officers as decision makers. The program has been well received by the policy makers in 
the USAF HEART program office and is being used as part of the long term project planning 
currently underway in the Air Force. For more details on the interactive computer program 
and its output, the reader is referred to Sadagopan [16]. 

6. CONCLUSIONS 

In this paper we have developed interactive procedures for the solution of bicriteria pro¬ 
grams. Even though the entire set of efficient solutions can be generated and pictured, such 
iteractive procedures do have merit due to the following reasons: in many real life problems 


f\(x) 

Mx) 



aj < xi < b, 

7-1, .... 17. 
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complete generation of all efficient solutions may be computationally expensive; even granting 
such generation, the DM is unaided in the selection process. Hence, the procedures developed 
in this paper may be thought of as "decision aids" [15]. 

The primary power behind our procedures is the unimodality property (Theorem 3). A 
similar result is available in Geoffrion [7], However, our proof of unimodality is entirely new 
and allows extension to more general cases unlike Geoffrion’s result. Also our relationship 
between Lagrange Multipliers and local tradeoffs is new. For more details on the bicriteria 
results, U.S. Air Force application and extensions to the general multi-criteria problem, the 
reader is referred to Sadagopan [16). 
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ABSTRACT 

There are a greai number of queueing x> stems. including the A/V.W 'i. the 
07 VA//1 and the discrete 07/0/1 queue in which ih te probabilities are 
determined h> repeated queue equations This pape- . .1 simple. eBicient 
and numerically stable algorithm to calculate the stir. ; ,i\ ' /ties and measure 
of performance for such systems. The method avt.i - in complex arith- 
mciric and mairis manipulations 


1. INTRODUCTION 

This paper gives a new method to find the equilibrium probabilities of continuous-time 
Markov processes in which the rate of going from state ttt to n is equal to d„- m for any n 
exceeding a certain limit, say r. The largest and smallest possible jumps are li and -g respec¬ 
tively. that is, the subscript of d k runs from -.g to h. These concepts will be further expanded 
in the next section. Here, we demonstrate it shortly, using the M x /M } /1 queue as an example. 
In this queue, arrivals of size ./', 1 < /' ^ It occur at rate A,. and departures of groups of size / 
occur at a rate of /a,. 1 ^ < g. Hence. 

r/, = A,. !<./</,■ 

(L, = n,. J < a. I < i < K 

The paper also considers discrete Markov processes with a similar structure. In this case, 
there is a probability of p k to go from slate n to n + k. Queues with such structure include the 
system size of the Gi x /M/c queue 14] and the waiting time of the discrete 07/G/1 queue I5j. 

At present, there are several approaches available to solve such queueing problems. In 
the first approach, one sets up generating functions, and uses partial fraction expansions to 
invert these generating functions. These expansions require one to find the roots of an equa¬ 
tion of degree g + h. The same is true if one uses operators on the queueing equations and 
determines the characteristic roots. I or practical applications, these approaches have some 
disadvantages. The equations are often difficult to solve, especially when they are of a high 
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degree. Their roots are usually complex and difficult to work with. Moreover, there mat be 
double roots which lead to further complications Neuls (3l therefore suggested solving these 
problems using matrix equations. 

The method in this paper combines the approach of Neuts with the classical method. In a 
sense, it can be considered as a modification of the method of Neuts. The main difference 
between our method and the Method of Neuts is the following. In Neuts' approach, a matrix R 
of dimension h x /; has to be calculated. In our approach, only the first column of R is calcu¬ 
lated. which reduces the number of operations approximately by a factor of h. Even though the 
method is similar to the one of Neuts, its derivation is done, using the classical approach. In 
this way, our method bridges the gap between the classical approach and the approach of Neuts. 


2. PROBLEM DEFINITION 

In this section, a precise problem definition is given and demonstrated, using a number of 
examples. First, we discuss the discrete case. Let there be a Markov chain with the transition 
probabilities. n.m ^ 0. The equilibrium probabilities of this chain are denoted by ir„. 
that is, 77 „ is probability to be in state n. As is well known, the ir„ arc determined by the fol¬ 
lowing equilibrium equations, provided the system has an equilibrium 

(1) 77 „ = £ Jr m />„,„• n > 0 . 

We assume that there is a certain limit rsuch that for n ^ r. p m „ = only depends on the 
difference between m and n. In other words, if Q represents the state. d k is the probability to 
increase Q by k from n to n + k. If ir„, is defined to be mo for m ^ 0, equation U) becomes 
therefore for n ^ r 

77 „ =» £ Vmdn-m- £ 77 „^* (i k . 11 > 7. 

We furthermore assume that d k = 0 for k < -gor k > h. Thus, one has 
I » > r. 

If one solves this equation for 77,,, one gets 

(2) ir„ = £ n„.. k e k . n > r. 

The e k in this expression are defined as 

e k => d k /(\ - d t) ). **0] 

(3) ,o-0. | 

For n < r. one uses equation (1), that is 
14 ) 77 „ — £ 77 mPm.n- n < r - 

Equation (2) will be referred to as the queueing equation , whereas equation (4) gives the initial 
conditions. There are r initial conditions. However, it can be shown that one of them is redun¬ 
dant. In its place, one has the normalizing condition 
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( 5 ) 

IT there exists a unique equilibrium vector [w„l. this vector is fully determined by the queueing 
equation, the r — 1 initial condition and the normalizing condition. Of course, if r = 1. there 
are no initial conditions. As is well known, an equilibrium solution will exist if 

p- iM/i< i. 

This paper will show how to find this equilibrium solution in an efficient way. 

In continuous-time Markov chains, the equilibrium probabilities are ueiermined by 
<6) 0 - £ it m q m „. n > 0. 

Here. q m „. m & n. are the transition rates, and 

We now consider problems where 

q m „ - d„- m . n > r. 

Here. d n is defined as 

rfo-- I<4- 

As in the discrete case, it is assumed that d k = 0 outside the range —y ^ f < It. If n„ = 0 for 
n < 0, one finds thus from equation (6) 

0= £ ir„- k d k . n > r. 

This equation can again be solved for rr„, which gives 
(7) £ *,- k e k . r. 

Here, 

e„ = 0 J 

<8) t\ = d k /l-d„ ). fc ^ O.J 

Note that equation (2) and equation (7) have exactly the same structure. In particular. e k ^ 0 
and the sum of the e k is one in either case. The normalizing condition is given by (51 both in 
the discrete and the continuous case. The initial conditions of the continuous case are easily 
found from (6) by letting 0 < n < r. 

We now give 3 examples which will be used later for the purpose of demonstration. 

EXAMPLE 1: First, consider the waiting time in queue of a discrete G//G/1 queue. Fol¬ 
lowing Ponstein [5], a Markov chain is set up, in which the state space are the waiting times 
rather than the number of elements in the system. To set up this chain, it is assumed that the 
difference between the service lime (5) and the interarrival time <>1) is (a) always integer and 
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tb> always between -g and li. It is well known (see e.g. Ponstein 15]). that the wailing lime of 
customer / exceed the one of customer / — I by an amount .S' — .4. provided, of course, that 
customer / - 1 had a waiting time above — (S - 4 ). Consequently, if W represents the wait¬ 
ing lime of a customer in equilibrium, W increases by A with probability P(S — ,4 = A). This 
means 

d K = P(S - A = A), -g A ^ h 
and 

e k = P(S - A = A)/[ 1 - m = A )]. A * Ol 

1,1 *-• I 

In the problem considered, there is only one initial condition, and this condition can be 
dropped as redundant. The queueing equation, together with normalizing condition is thus 
sufficient to determine all ir„. Here, 7r„ represents the probability of having a waiting time of 


EXAMPLE 2: In the Gl h /M/\ queue, arrivals occur in groups of size h. Let a , be the 
probability of servicing / elements between two consecutive arrivals and let rr„ the probability 
of having n elements immediately before an arrival. It is easily verified that the queueing equa¬ 
tion is given by (2) with r = 1 and 


e k = a h .. k /d - a h ), < A < /;. A ^ 0| 

^ ,o=0. ( 

Here. a k is the probability of serving A customers between two successive arrivals. The initial 
conditions consist only of one equation which can be omitted. Also, in the case considered, g 
is infinite. However, for practical calculations, it is sufficient to carry only a finite number of 


EXAMPLE 3: Consider the M x / M y /1 queue. A ( , 1 < ./ < /?. is the rate at which groups 
of size ./ arrive. No service is done unless at least /.’elements are in the system. If the number 
in the system is g or more, there is a rate of (i,. I < / < g , that i elements leave the system 
after having received service together. For this problem, one has the fQllowing transition rates 

(/„.„ = \„- m . 0 < // — m $ h 

dm.,, = Pm-n- 0 < m - « ^ g, m > g. 

All other q m n . m ^ n are zero, in particular, all m < g, n < m are zero because nobody 
leaves if there are less than g elements in the system. The q m m become 

q„ m = - 21 A , m < q 
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A little reflection shows that r = g. and the queueing equation is given by (5) with 


e„ = 0 



| £ x, + £ m, 

|. A > 0 

(11) 

G = Ma/ 

[|>/ + i>j 

. A- < 0 

The initial conditions are given by (7), that is 


0 = £ » = -7r <> £ K > + 

o = £ = £ ~ n " £ A ' + £ « < & 

Again, the convention was used that n m = 0 if m < 0. From the initial conditions, one can be 
omitted, and it turns out to be convenient to omit the one corresponding to n = g - 1. 

Above, we discussed a few examples which are amenable to our method. Additional 
examples can easily be generated. We now show how to solve these examples. 

3. THE MAIN RESULT 

The fundamental result of this paper is the following. If p < 1. and if the Markov chain 
is ergodic, there are h values a, > 0, 1 ^ j < h such that 

( 12 ) n n = £ a ; 77„_„ n > r. 


The a r 1 < j ^ h are determined in a unique way by the e k . -g < A < /;. Further¬ 
more, 

£ a, < 1 . 

This paper will present several methods to find the a,. Once the a, are found, the problem is 
actually solved. Together with the r — 1 initial conditions and the normalizing equation (12) 

gives exactly r + 1 equation io determine rr 0 . ff|. n r . Moreover, the generating function 

of the 7 t„ can be expressed in terms of a and no. n t .To see this it is convenient 

to define a n as -1, and set 

A (z) = £ a.z' 

/*<*>- £ rr„r". 

in this case, one finds, as will be shown in Section 4: 

(13) /l(z>P<z) = £ o,£ w^z*. 
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When setting z = 1. £(z) becomes I, and 

- 4(1)= £ a > I £«,. <' = min(r - n - 1./;). 

This means 

(14) £ a, = -1 + a| + a 2 + • • ■ + a h — £ ir, £ a,, e = min(/- - « - 1, 

This equation provides a convenient normalizing condition and can be used in place of (5). To 
find L. the expected number in the system, one can differentiate (13) and find 

A'(z)Piz) + Aiz)P'(z) - £ a, £ nn 

If z — 1, P(z) = 1 and £’(z) = L. Using these values, the above equation gives after some 
calculation 

"51 L- 3>,|l jj / j>. 

Equations (14) and (15) simplify if r = 1. Then 
A (z) P(z) = = —■wo¬ 

lf one sets z = 1, this gives. 

(16) ir n = -A(l) - 1 - £ a,. 

This means that 

I a, < 1. 

If Q is a random variable representing the state of the system, £(£)) can be found as: 

(17) £((?) = £ nir„ = -A'(l)/A(\) = £ ja,/n () . 

The ideas just presented shall now be demonstrated, using the examples discussed above. 
First, consider the discrete Gl/G/l queue. Specifically, suppose that the distribution of interar¬ 
rival time and the distribution of the service time are as given in Table I. 

T ABLE 1 — Distribution of the Interarrival Time 
and the Service Time. 



Using these values, one can find the e k , -g ^ k < h, from equation (9), and the e k in turn, 
determine the a,, 1 < j < h (see Sections 4 and 5). For the problem described by Table 1. 
one finds, 

a, = 0.23156 
a 2 = 0.05039. 
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Equations (16) and (17) give now 

- PiW- 0) = 1 — a, — fl2 = 0.71804 
E ( W) - (a, + 2a 2 )/P( W = 0) = 0.46287. 

Since w.) is assumed to be zero, one can use (12) to calculate 7T|, n 2 , jr ; and so on. 
ir i = 7T ()& | + 7r_|fl 2 = 0.1663 
n 2 = 7T i zz [ + 7T|)ai = 0.0747 


Since all a, > 0, this recursion is numerically very stable. 

Next, consider the Gl h /M /1 queue. Specifically, let the arrival time distribution be deter¬ 
ministic, let h == 2, X 2 = 0.4 and n = 1. The a k of equation (10) are in this case Poisson- 
distributed with parameter 2.5. One finds in this case 

o, - 0.41045 
o 2 = 0.13711. 

Again, all entities of interest can be calculated. L. the expected number in the system preceed- 
ing an arrival becomes for instance 

L - (oi + 2flj)/(l - a x - a 2 ) = 1.5132. 

Finally, consider the M x /M y /\ queue of Example 3. Specifically let 
= \ 2 = Mi = M2 = Mi = *• 

In this case, the aj turn out to be 
o, = 0.34960 
a 2 - 0.24582. 

The initial conditions are now 
0 - -2n 0 + ir., 

0 = 7r 0 - 2n i + ttj + it 4 , 
ir } and ir 4 are given by ( 12 ) as 

iTy = 2 + « 2 n 'i = 0.349607r 2 + 0.24582ir t 

tt 4 = a t 7 r 3 + a 2 n 2 = (a j 2 + a 2 )ir 2 + a,a 2 ff| — 0.368047 t 2 + 0.08594^). 

Using these values, the initial conditions become 

0 - -2tt 0 + 0.24582tt, + 0.34960ir 2 
0 — 7r n — 1.66824ir,+ 0.71764 tt 2 . 

As a third equation, one uses the normalizing condition as given by (14). 

— 1 + a, - 1 - a 2 — ito (—1 + a i + a 2 ) + ff| (—1 + a t ) — n 2 
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0.40458 = 0.40458tt o + 0.65040,7, + ir 2 . 

These 3 equations ean be solved for 7 r () . ir\ and 7 r 2 . 

The result is 

7 T„= 0.06741. 7r| = 0.15840, n : = 0.27428. 

It is now simple to find L from Equation (15). One has 

L = I— <7r( + 27 t 2 ) + a 1 (jr,) + 2 tT| — 1) + a2 (2n n — 2)1/(— 1 + a, + a : ) 
= 3.4128. 


Thus, once the a, are found, the problem is solved easily. The next two sections show 
several methods to find the a,. 

4. THE SIGNIFICANCE OF a, AND THEIR CALCULATION 

In this section, it is shown that equation (12) holds, and how one can find the a, in the 
general case. From the initial conditions and from the queueing equation, one finds the follow¬ 
ing expression: 

Piz) = Uiz)/Viz). 

Here Uiz) is a polynomial of degree r + g - 1 which is of no further interest, and l'(:l is 
equal to 



It is well known (see e.g. |5l) that U(r) has a zero at r - 1, j? — 1 zeros inside the unit circle 
and h zeros outside. It follows that V(z) can be written as 
Viz ) = Aiz) Biz). 

Here 

Aiz) — —1 + 0]Z + a 2 2 2 + ... + a^z 1 ' 

contains all the h zeros outside the unit circle. 

Because of the convergence of Piz), Uiz) must be divisible by Biz). Since Uiz) has a 
degree of r + a - I. and Biz) has a degree of n. Ciz) has a degree of r - 1. Consequently. 

Uiz)/Biz) = Ciz) = £ c,z‘. 


Piz) = Ciz)/Aiz). 
This equation can be written as 
(18) Piz)A(z) - Ciz). 
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The coefficient of of the product P(r).4 (z) equals, as is well known 

'L a '”»-r 

for n ^ r - 1. this means 

(19) £(t,7T„_,= C„. 

For n > r. one finds 

(20) £a ( ir„., = 0. 

Equation (20) is identical with Equation (12), and Equation (18) and (19) can be combined to 
give Equation (13). This proves Equations (12) and (13). 

From the above discussion, the significance of the a, becomes clear. I (c) can be factor¬ 
ized into two factors .4(c) and Biz). and the a, are merely the coefficients of 4<c>. provided 
is defined as -1. The problem is thus reduced to finding 4(c). The most convenient 
method to do this seems to be the following one. One starts with the identity, which will be 
proven in Section 6. 

(21) a, = e, + £e_, «,../• 1 < ./' < /'• 
with 

a () , - a, , I <./'</) 

<22) a,*!., - o, | a, + a, l+l , I ^ < g - 1. 1 j < h 

Furthermore. <?,, is to be taken as zero fo. j > h. Equations (21) and (22) can be used to find 

a, by successive approximation. As starting values, one can use a" = e,, j — I, 2. Ii. 

These values can be used to replace the a, on the right of (22), and a new approximation for a, 
can be found from Equation (21). In this way, one continues until a suitable stopping criterion 
is satisfied. 

We tried two different stopping criteria. First, we required that the it, generated by Equa¬ 
tion (12) satisfy the queueing equation with a specified precision e. It was found that for high 
values of h and/or & this criterion performed poorly. The reason for this is simple. If h is say 
100, a small change in d m will have a dramatic impact on p and, consequently, on the perfor¬ 
mance measures of the system. As an alternate stopping criterion, we used the change of the 
calculated .4'(1) = I ja, to decide when the precision of the a, is adequate. Since 4(1) is 
closely related to the expected number in the system (see Equation (17)) this seems to be a 
good stopping criterion. 

If there are q iterations, one needs 4qnh operations to find all a,. The reader may want to 
verify this. If one adds all e_,a, , to the sum at the right of Equation (21) as soon as the a, , 
becomes available one only needs to store the a u for the current value of i. If this is done, the 
algorithm requires only an array area of 4h + #as the reader may verify. The algorithm is thus 
efficient and does not require huge arrays. 

The algorithm described above was programmed in FORTRAN and run on the DEC 
2060. The execution times were negligible. A problem with h - 100 and g - 100 and p = 0.9 
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«»nK required X 4 seconds and 53 iterations until all a , were obtained. At this point there was a 
relative change of A '< 11 of less than 0.0001 between two iterations. The computer time to find 
ihe a lor e - 100 and ft = 5 is 0.2307 seconds (19 iterations). The method performs thus 
vetv well 

5 ALTERNATE METHODS 

In gam a check for our results, and also to get further insight into the problem, we 
decided h> try other methods as well. These methods are based on IT;) = L(z)/(z - 1), that 
/is I i;> was deflated by ; — I. IT;) can be written as 
It;) = 4(;)fl(;), 

where H(:> is equal to #(;>/(; — 1). This deflation decreases the degree of the polynomials 
one works with by one. and it also increases the difference between the zeros of T(z) and 
H (;i. improving thus the efficiency of the algorithms to be discussed. A further deflation by 
the only positive zero of A (;> is possible. 

To factorize It;), two methods were used. The first one is given in U, page 158). The 
second one used the fact that IT;) can be interpreted as the characteristic polynomial of a 
difference equation If these difference equations are used to calculate values x„, x / ,+ \, jc „+2 
recursively, the zeros greater than one will dominate, and eventually, the effect of B(z) will 
become negligible The ,v„ are thus almost identical to a series generated by a difference equa¬ 
tion that has ,4(z) as characteristic polynomial. In other words, the x„ will almost satisfy the 
following relationship 

If x„ is known for n — A, k + 1. k + 2ft — 1, this gives h equations for the h unknown 

a r 

The x„ will satisfy Equation (23) precisely if the initial conditions xq. xi.x*_| satisfy 

(23). Consequently, one can repeat the above algorithm several times to gain a higher preci¬ 
sion for the a r In each repetition, one uses the a, obtained in the previous iteration in order 
to calculate the starting values needed. In the first iteration, one can use a,° - d r This algo¬ 
rithm gave good results for most cases we tried. 

The disadvantage of the two algorithms just mentioned is that they require the solution of 
h equations in h unknowns during each iteration, and this requires h } operations. This means 
that they are inefficient as compared to the algorithm suggested in the previous section, at least 
if h is high. 

Of course, A (;) can also be found, provided one knows all the zeros of T(;) outside the 

unit circle. Let Z|,Z2.zj. z h be these zeros. A(z) becomes, given one uses the fact that 

a o - -1 

A (;) = -(; - z\) (z - z 2 ) ... (z - z*)/[(-Z|) ( -z 2 ) ... -(-z*)] 

= -0 - z/z\) (1 - z/z 2 ) ... (1 - z/z h ). 

The factors can now be multiplied in the usual way, (see e.g. (4J) giving 
A (r) — — 1 + a\ z + fljz 2 + ... + a h z h . 
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Note that this procedure even works for multiple roots, provided, of course, all are known. 


A final comment shall now be made. Suppose one has the h + g probabilities 7r**|. 

ir A+: .7 T k+h t x - Then, one can theoretically calculate ir m , m > k + h + g. recursively by 

using Equation (2), that is 




This recursion is numerically unstable, that is, the round-olT errors will increase with each itera¬ 
tion. The characteristic equation of this difference equation is fd/.v) = ,4(l/.v> Stl/.v), and 
any recursion based on it will eventually be dominated by fl(l/.v). (Indeed, this very effect 
was used earlier to find /Mr)' Instead of doing such a recursion, it is better to use 2h subse¬ 
quent 7r„, and use (12) to obtain //equations for the a,, j = 1.2.//. 


6. THE SIGNIFICANCE OF THE a u 

This section proves Equations (21) and (22) and establishes the relationship of our 
method with the method of Neuts. We start to prove the following equation, in which the a ,, 
are calculated as given by (22) 



Since a {) , = this equation is certainly correct for / - 0. Moreover, if it is true for i. it is 
also true for /-El. One has, if a, h+i = 0 

7r, l4 .,4.| = £ a,, 7r, I+ i- ; = a lA ir„ + a )( jr„+i-, 

= «,, £a,7r n .,+ £fl, )+l rr„_, 

= £ la, |U, + a,. /+ |]7r„_ ; . 

Because of (22). this gives 

*V.+ ! = £ a i+\.i 1 r n-r 

This proves that Equation (24) is correct for (+1, given it holds for / To prove Equation 
(21), one rewrites the queueing Equation (2) as 



Using Equations (12) and (24) this gives 


I 


£ + £ e -> Z a ' i ir n- l 

£ [e, + c-, 
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Consequently, 

a i = e i + Z e ~< a t r 
and Equation (21) is proven. 

Equation (24) can also be used to simplify the initial conditions. Consider for instance 
Equation (4) which can be written as 

lr„ = £ TtmPm.n = £ W n,Pm.« + Z Pr+m. n "r+ m - » < r - 

If n r+m is replaced by (24), one obtains after some calculation 

ir„ = £ rr„ + £ n,„lp m „ + £ n r *k.„ a kr _ m ], 0 < n < r. 

Together with Equation (14), this provides a set of r independent equations in the r unknowns 

77(|, 77 l.77,_). 

We now compare our results with the ones obtained by Neuts. [3,41 Neuts sets up the 
following matrix equation for the unknown matrix R. 

(25) R=^R"A n . 

The A„ are problem-dependent matrices. In our case, they become 

■4,,= [4*1. if-0.1.2. [gfh 1 + i 

Here, [#//;] is the lowest integer above j?//;, and 

A,", = 

All other A„ are zero. As before, all e k outside the range -g < k < h are assumed to be zero. 
If one defines 

77* = [77 A „ + |, 77*7,+ 2 .7r*/, + /,i. 

one has according to Neuts 

77* + l = 77*/?. 

If R = [r, ; ), this gives 

or 

w B+( - 21 7r «+i- f O.+i-,.i- 

When this equation is compared with (24), one finds 
In particular. 


a 0.i “ Oj — r h + 1-/. I- 
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The a , are thus contained in the first column of R. We also note that Equation (21) follows 
from Equation (25). as the reader may verify. 

To calculate R according to (25) one has to find the terms R"A„. n = 2.3. \nlh\ + I 

and this can be done in 2U///1 matrix multiplications. However, each matrix multiplication 
requires 2h‘ operations, which means that there are 4 /r' [#//;). or approximately 4 h : ,v opera¬ 
tions per iteration to find R . not counting the matrix addition Thus, by taking advantage of 
the special structure of the problem one obtains a more efficient algorithm, namely the one 
described in this paper. 

The equivalence of the a, with r h+ 1., | is of great importance. In particular. Neuts (21 
proved that all r u are nonnegative. Consequently, all a , are nonnegative. Neuts also finds that 
the r,, have a probabilistic interpretation. Indeed, the r, , are the expected number of visits to 

state n + H + i starting at n + I without hitting (// + I. n + i . n + //). This means 

that a, gives the expected number of visits to slate n + / starting at n without hitting 
I n + ./-/). n + ./ - 1). 

7. CONCLUSIONS 

This paper derives an extremely efficient algorithm which can be used to numerically 
determine the slate probabilities of many queueing problems efficiently and precisely. This 
algorithm was obtained a by combining the classical approach with the approach of Neuts. 
Indeed it bridges the gap between these two approaches and opens new horizons for further 
research. 
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ABSTRACT 

This paper examines the process by which a user of a queueing system 
selects his arrival lime to the system to compensate for unpredictable delays in 
the system if he wishes to complete service at a particular time. Considering 
the case in which all the system users have already decided on their arrival 
times to the system and will not change these limes, this paper investigates how 
a new user of this system develops his strategy for selecting his arrival time 
The distribution of this customer's arrival lime is then obtained for a special 


I. INTRODUCTION 

If a user of a queueing system can choose when to enter the system, then he would prob¬ 
ably choose to enter when the delays he experiences will be minimum. However, if he wishes 
to complete receiving service at a particular time, called the "target time," then he would have 
to arrive at the system early enough to allow for his service time and unpredictable delays in 
the system. This allowance was termed "headstart" by Gaver (3) and the same terminology will 
be adopted here. The strategy by which a user selects his headstart is studied in this paper. 

Consider the case in which all the system users already have fixed headstart strategies and 
will not change these strategies under any condition. This paper investigates the case in which 
a new user who wishes to use the system develops a headstart strategy. This will be done by 
assuming that this user, whom we name U (cf. Alfa and Minh [2D, attaches some perceived 
costs to delays and to early and late completion of service, and that he selects his headstart in 
order to minimize his total cost. Gaver [3] obtained the "best" headstart strategy for this user, 
in this type of situation, on the assumption that the user wishes to minimize his expected total 
cost. This paper proposes the headstart strategy usually adopted by such a user, on the assump¬ 
tion that he selects his headstart each day in order to minimize the total cost he incurs. It is 
assumed that on the first day he uses the system he does not know what the delay distribution 
is and therefore selects his headstart arbitrarily, and hence probably incurs high total cost. 
However, for the following day, he uses his knowledge of the previous day's outcome to try 
and reduce his total cost; and this continues until a stage is reached in the long-run when he 
settles for a particular headstart and does not change any further. Let us call this stage 
"steady-state." This headstart reflects his arrival time at the system and our interest is in the 
distribution of his arrival time at the system at steady-state. 
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This paper does noi restrict its analysis to nonrush hour traffic or linear cost parameters as 
done by Gavcr [3]. 

2. PROBLEM FORMULATION 

2.1 Assumptions and Definitions 

Consider a queueing system used by a finite number /. (/ > 0), of customers who behave 
identically and independently. For each day. let us observe the system at equally spaced epochs 
sequentially numbered 0.1,2,3,... and assume that all the customers arrive at and depart from 
the system only at instants immediately prior to these epochs. 

It will be assumed that once a customer is served he does not return to the pool of poten¬ 
tial customers. This assumption and that about the number of customers being finite are not 
necessary for the analysis but they allow us to use the results from Alfa (I] and Minh [4). 

Assume that each customer, considered separately, has the time-dependent probability K„ 
of arriving at the system prior to epoch n + 1, given that he has not arrived at the system by 
epoch n (// = 0,1,2.3,...). Assume that the service times of the successive customers are 
independent, identically distributed random variable S in the set {1.2. MY. M < 

Suppose the new customer. U who also wishes to use this system, chooses to arrive at the 
system at epoch a = n. Let the delay he experiences by so doing be W". The distribution of 
W" can be obtained as shown in Alfa [1). Note that the addition of this customer U to the sys¬ 
tem has raised the total number of customers to / + I. Although the distribution of l/s arrival 
time is not necessarily the same as that of the other / customers, the results in Alfa [1] can still 
be applied keeping in mind that W" is a conditional random variable. The distribution of I/'s 
arrival time shall be developed in this paper. 

Let V" be the time spent in the system by U, given that he entered the system at epoch «, 
then V" = W" + S. Let V" A Pr(V” = /). 

For practical purposes it will be assumed that U shall arrive at the system only at an epoch 
between 1 and ,V; where I < S < °o. 

Suppose U wishes to complete service at epoch r, where for convenience 1 < r < ,V. If 
he arrived at the system at epoch a = «, got delayed W" = / units of lime and took S = / units 
of time for service then he would complete service at epoch n + i + j. One of his objectives is 
for n + i + j to equal t. However, if n + /'+./< r or n + /+./> r then he is early by an 
amount r - rt - / - / or he is late by an amount « + /+./'- t. respectively. He attaches a 
cost to either of the outcomes and these costs are not necessarily the same. In addition, he 
attaches a cost to the time spent in the system i + j. He would, therefore, incur a total cost 
C'(n. i + ./'), where C(n. i + j) is the sum of the cost of time spent in the system and the cost 
for either earliness or lateness, depending on which is the outcome, keeping in mind that if 
n + i + ./' = r then he only incurs the cost for his time in the system. 

Let C„ be the total cost incurred by U given that he arrived at the system at epoch n 

where C„ can only assume the values in the set - 1.2.(/ + II x M |. For 

brevity, C„ will be termed the cost associated with epoch n. 
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Let us define the following sets that will be used in the next subsection. 


(i) G A |tit = 1,2.1/ + I) x \1 1 

(iil L A Ii 11 = 1.2.3.V} 

Hit) L„ A (/I/ = 1.2.3. n - I. n + 1..V| 

<iv> LI A L„ n l.„. 

If there are v epochs. (1 < v ^ .V - 2), (.V > 2). labelled m t . m 2 .m, for which 

C m = C m| = C m , = ... = where Vm r € Z.„. r = 1,2.v. then let 

<v> L” ms A {/It - m : .« v ; Vm, € L". r- 1.2.v|. 

Let 

(vi) I", it, - L”„ i.e. L " v is the set difference of Z. m and JL£,. 

2.2 The Model 

2. 2.1 The General Model 

Suppose on the d" 1 day, (rf- 1,2, ...). Li arrived at the system at epoch a 1 '= n. Let 
alf A Pr(a‘ /= n|, n € L. U would incur a total cost C„. and by definition 

(1) Pr{C„ = C(n. /)) - V”. 

We assume that ITs choice of arrival time epoch for one day depends entirely on the out¬ 
come of the previous day's choice of arrival time i.e., on the total cost incurred. U wishes to 
minimize his total cost. Hence, a‘ <+l is a decision variable with an action space L. and a* 7 * 1 
depends on the cost incurred on day d. Given the action on the (f h day the action for the 
(d + l) ,h day, a d+1 , is thus a Markov Chain. 

Let t„ m A Pr(a J+l = m la‘ , = n) be the transition probabilities of the chain. These transi¬ 
tion probabilities define If s strategy. 

U s intention is to choose an arrival time epoch that minimizes his total cost each day he 
uses the system, hence for the (d + l) th day, his strategy for choosing his arrival time epoch 
will be 

(2) - Pr{C„ < C„; V« 6 L „); 

(3) i„. m = Pr{C m < C H ; Vw € Lj 

+ I - Pr|C„, = C„ r < C„; Vm, € Li.,. : Vw € L", J; 

Vm n . 
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The arguments leading to (2> and (3) above are as follows: 

(a) If L arrived at the system at epoch n on the tf h day he will not change his arrival 

time to any other epoch on the id + I ) lh day if he can not reduce his total cost by so 

doing. This leads to (2). 

(b) If U arrived at the system at epoch n on the </"' day he will change his arrival lime to 

epoch m on the id + l)' h day if this could reduce his total cost, provided there is no 
other epoch which has equal potential reduction in total cost as epoch m— this leads 
to the first term on the R.H.S of (3). In addition, however, if there are v other 
epochs, other than n and m. at which L could choose to arrive at the system and 
achieve the same reduction in to’.;.! cost as arriving at epoch m. then L would arrive 
at any of these v epochs or epoch m with the same probability, provided the total 
costs to be incurred at each of these epochs, including epoch m. is less than at any 
other epoch. This leads to the second term on the R.II.S. of (3). <2) and (3) above 
constitute the transition probabilities that define L's strategy. 

Define an ,V x V matrix T = (/„„,) and an ,V vector \ J = Ini', a-. ai I. Then L 

selects his id + I) lh day's headstart according to 

(4) A 1 '* 1 = A 1 ' x T. 


THEOREM 1: £ i„ m = 1. 

PROOF: Consider one epoch n, in € L), and the associated cost C„. If we compare this 
cost C„ to the costs associated with other epochs in the set L„ then C„ is either the minimum 
cost, one of the equally minimum costs or neither of the two. This can be stated as 

(5) Pr{C„ C u ; V M € L „) + Pr(C fl <C„; V« € L m \ = I. 

Substituting i„„ for the first term on the L.H.S. of (5) gives 

(6) i nn + Pr(C„ <C„: Vm € L„) = 1. 

The second term on the L.H.S. of (6) is the probability of having at least one epoch, in 
the set whose associated cost is less than C„. This implies that there is at least one epoch, 
other than n, whose associated cost is minimum. If there is only one such epoch then the pro¬ 
bability of such event is given by 

£ Pr(C m < C» ;Vh’ € L m \. 

*m(l n 

If. however, there are v other epochs (I ^ v ^ .V - 2). labelled ni\. nt-_ .w v . Vm, € 

I < r < v. such that C,„ = C M| = C,„, = ... = C m> , then the probability of such event is given 
by 

v I I ttt Pr ( c - = c ’ m . < t »: v„, r e v» e r m t ). 
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The second term on the L.H.S. of (6) can thus be written as; 

(7) Pr{C„ <C„; Vu € L„) = £ Pr{C m -< C„; Vw € L m 1 

V-2 i 

+ £ £ y- Pr|C m = < C„; Vwt, 6 I». v , Vw 6 L ", v ) 

= v I 

Substituting this into (6) gives 

+ Sl = ' 

i.e., £ t nm = 1; and this proves Theorem 1. 

REMARK: If i„„ = 1 for any w € L then w is an absorbing state and the solution to 
Equation (4) is trivial. However, the occurrence of such a situation in real life would be quite 
rare —we therefore assume for our present problem that 0 ^ r„ m < 1, V(w,w/) € L. 

LEMMA 1: If t ni .„,= 0, <wi.w>) 6 L, for any n t ^ w 2 and r H| , (| ^ 1, then: 

(a) either t n> „ 2 = 0, Vw, € L, and w 2 is a transient state which can be deleted from the 
chain, or 

(b) there exists at least one w 4 w, such that i„ t > 0 and therefore state n 2 can be 
reached from state w \ via state w 4 , hence, all w € L form an irreducible markov chain 
(Note that w 4 6 L). 

PROOF. Lemma Ha) is a compliment of Remark above and docs not require a proof. 
For Lemma 1(b), if there exists w 4 ^ such that t n > 0. then there also exists at least 

one N tuple U],i 2 ./\) such that (~(w 4 ./ 4 ) > ('(wi./M ^ C'(/i,./,) ^ 

C(»i./|) < ('<«*,4); V* > 4. V/ ; e C, K j < N. 

Further let 

(vii) L 0 A {/14, = 0; Vw € /.} 

(viii) L 0 & L — L a 

<ix) L (l „ A L 0 n L„ and let yV 0 be the number of elements in L 0 . 

Define an A 0 x /V 0 matrix T 0 = U „ m ), V(w,wt) € Z. 0 i and also define an Ao vector 
A; - l afl: Vw € Z,„. Equation (4) can now be stated in a modified form, for only the epochs 
"rnained in L„. as; 

*■ A,'/* 1 = A,f x T 0 . 

I MM \ 2 i„„ > 0, Vw € E,,, hence all w € £ 0 are aperiodic states. 
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PROOF: From ihe cost structure, it is apparent that if t„ „ = 0 then r m „ - 0, Vm 6 L , 
and hence n ? L n i.e., state n would have been deleted from the chain considered in equation 

(8) (i.e., in L 0 ). Therefore, t„„ > 0, V« € L 0 , and hence Lemma 2 follows. 

THEOREM 2: There exists a limiting probability vector A 0 = lim A 0 rf such that 

(9) A 0 =A 0 xT 0 

has a positive solution which is unique. 

PROOF: The necessary and sufficient conditions for the existence of a positive and 
unique solution A () is that the chain should be irreducible and aperiodic. Both conditions were 
established by Lemma Kb) and Lemma 2 respectively. 

The interest of this paper is in LTa choice of arrival time in the long run, A 0 , which can 
now be solved for in Equation (9). 

2.2.2. A Special Case 

In the general model it was assumed that if U arrived at epoch n on the day, he would 
not change his arrival time to any other epoch if C„ < C m , Vm € L„. Suppose we modify this 
assumption and now assume that when C„ - C m for any m ^ n we would not rule out the pos¬ 
sibility of [l changing his arrival time to epoch m. In that case we shall attach equal chances of 
U changing his arriv;.; i'me to epoch m and to him not changing his arrival time from epoch n. 
Thus, If s strate 6 , will now be modified such that the transition probabilities that define his 
strategy will be given by; 

(3a) t n . m - Pr{C m < C„; Vu € L m \ + 

I Pr{C m = C„ r < C„; Vm r 6^; 

Vw 6 L nv ): Vn.m € L, 

where 

(x) L m v A [i\i — m\,m2 .w v ) such that 

c m - c m| = C m2 “ ... = C„ v ; 1 < V < N - 1. and 
Z. m . v c L \and 

(xi) L m v A L m - L„ v . 


It is immediately apparent that the right hand side of Equation (3a) is independent of n, 
hence i„ m = r vm ; V(n,m,\) € L. Let t m A t„ m \ \{n.m) € L. For this special case, therefore, 
all the rows of T are identical, and the steady state solution is given by 


( 10 ) 


/„ ; Vn 6 L„ 
a " = 0 V« 6 L„. 
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2.2.3. Illustrative Example 

For illustrative purpose let, 

/ = 2, N- 3, r = 3, M = 1. 
\„= 0.3, A, = 0.5, \ 2 = 1.0. 


Let; 


( 11 ) 

and let C\. 


C'(n.i) = C u 


</) + 


(r - n - i) . 
(// + i - t) 


1.0. C, = 1.0, C„ = 2.0. 


By using results from Minh [4] we obtain the distribution of queue lengths and from Alfa 
[1], the distribution of waiting is obtained as 


W’’ 



The transition probabilities are thus given as 

(12) t ,, = f , 1 [ v\ + v\ ] + v\ [ v\ + j [ f 2 3 + vl ] + f , 1 v} vj. 

(13) / 12 = F, 2 + F 2 2 F, 1 [ F 2 3 -f F 3 1, 

(14) r ! 3 = F, 3 [F 2 ‘ + F, 1 ] [F 2 2 + F, 2 ] + V} FjF 3 2 , 

(15) / 2 ,i = Vf tF 2 2 + Fi) IF 2 3 + Fj 1 ) + Fj [F 2 2 + Fj 2 ] [F 3 + Fj 1 ] 

+ Ki 1 F/ Fj 1 + F,'[F 2 2 + F/]F|V2, 

(16) r 2 2 - F, 2 + F 2 Fj [ F 2 + Fj 3 ]. 

(17) r 2 , = F 3 Fj [F, 2 + F, 2 ] + [F, 3 + F 2 3 ] V} Fj + F, 3 F,' [F 2 2 + F, 2 ]/2. 
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(18) f 3 , = [ v} + v\ ][r? + v]\vl + [vl + v}][vl + v}Wl 

+ V}V}Yl 

(19) l J2 - V 2 + K 2 2 V\ [ Vl + Vj\. 

(20) = v} [u 3 2 + Vjt ) + I/, 2 ] + nU'/rj. 


Thus, this leads to the transition matrix, 
0.333 0.605 0.062 

(21) T = 0.256 0.605 0.139 

0.167 0.605 0.228 

from which we obtain 


(22) A-[0.265. 0.605. 0.130). 

Note that for this example L 0 - L, hence T 0 = T and A 0 


A 


3. CONCLUSION 


This model, based on its assumptions, can only be used to predict the customer’s arrival 
time at the system when all other customers previously using the system have settled for a par¬ 
ticular arrival time and will not change it under any circumstances. However, when all the 
other customers change their arrival times as a result of the new customer’s "interference" then 
there will be a slight modification to the problem. 

If other customers change their arrival times everyday, then the distribution of W" will 
change from day to day and, hence, T and T 0 will also change from day to day. If we let T rf or 
To' represent the transition matrix for the strategy for the (d + l) lh day then If s choice of 
arrival time can be reformulated as 

(23) A rf+I = \ d x 

However, with (23) existence of a steady state solution cannot be assured. If we further let 
there be more than one server in parallel such that a customer develops a strategy not just for 
selecting his arrival time only but for joint selection of arrival time and of the server to serve 
him, then the problem becomes similar to that in Alfa [2j. 
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ABSTRACT 

This paper develops a forward algorithm and planning horizon procedures 
for an important machine replacement model where it is assumed that the tech¬ 
nological environment is improving over time and that the machine-in-use can 
be replaced by any or the several different kinds of machines available at that 
lime. The set of replacement alternatives may include (i) new machines with 
different types of technologies such as labor- and capital- intensive, (ii) used 
machines, tiii) repairs and/or improvements which affect the performance 
characteristics of the existing machine, and so forth. 

The forward dynamic programming algorithm in the paper can be used to 
solve a finite horizon problem. The planning horizon results give a procedure 
to identify the forecast horizon T such that the optimal replacement decision 
for the first machine based on the forecast of machine technology until period 
T remains optimal for any problem with horizon longer than T and. for that 
matter, for the infinite horizon problem. A now chart and a numerical exam¬ 
ple have been included to illustrate the algorithm. 


1. INTRODUCTION 

In our previous paper 151, we developed forward algorithm and planning horizon pro¬ 
cedures for a machine replacement model under the assumption that only one kind of machine 
is available for replacement in any given period. We showed that there exists a forecast horizon 
T , such that the optimal replacement decision for the first period (based on the forecast of 
machine technology until period T remains optimal for any longer (than T) horizon and. for 
that matter, the infinite horizon. 

In this paper, we relax the assumption of a single possible replacement alternative by mul¬ 
tiple alternatives. That is, we develop a model in which the machine-in-use can be replaced by 
any of the several different kinds of machines available at that time. This model can deal with 
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many realistic situations. For example, the set of replacement alternatives may include III new 
machines with different types of technologies such as labor- and capital- intensive, (ii) used 
machines, (iii) repairs and/or improvements which affect the performance characteristics of the 
existing machine, and so forth. It should be noted that the set of replacement alternatives can 
change from period to period. 

This multiple alternative replacement model represents a major generalization of the vari¬ 
ous single alternative replacement models developed by Terborgh [71, Thompson [81, Sethi and 
Morton [61, Gordon 111. For a survey of the extensive literature on the machine replacement 
models, the reader is directed to Rapp [4); see also Pierskalla and Voelker [31. 

In the next section, we formulate the multiple alternative replacement model under 
improving technologies. In section 3, we develop an efficient forward algorithm for solving 
finite horizon problems. Section 4 develops a variant of the regeneration- monotoniciiy pro¬ 
perty [2, 5] and planning horizon procedures for the model. A flow chart summarizes the com¬ 
plete solution procedure. A numerical example is solved in Section 5 to illustrate the steps of 
the forward algorithm and the application of the planning horizon procedure. We weaken the 
assumption of improving technologies in Section 6 and show that planning horizon procedures 
can be adapted to this case. Section 7 concludes the paper with some important remarks. 

2. MODEL FORMULATION 

Consider the situation of a production shop which must keep a single machine of a partic¬ 
ular capacity at all times To run this machine, the shop incurs operating expenses which may 
include labor cost, electricity cost, maintenance cost, depreciation, and so forth. Usually, the 
performance of the machine deteriorates, i.e., operating expenses increase over lime, so that 
the shop might consider selling the existing machine for its salvage value and buying a new 
one. This new machine is selected from various different alternative machines available. If a 
major repair is decided on the existing machine, then this alternative can be considered as if a 
new machine of another technology is purchased at the cost of repair plus the salvage value of 
the existing machine. Once the new machine is purchased, the same situation repeats with the 
new machine. Consequently, for any finite or infinite horizon, the shop will make a chain of 
replacement decisions [4], [51. The problem of the shop, of course, is to find simultaneously 
the optimal times of these replacements and the types of machines selected at these times. 
Obviously, these decisions will depend on the prices of the future machines and their per 
period operating expenses over time. More precisely, we are considering the following model 
which the shop must solve: 

Let A* denote the machine of technology h available at time t. Let there be machines of 
N alternative technologies available to choose from in any period. It is, noted that the extension 
to the case, when the set of alternative technologies is changing over time, is straightforward 
and will be described in Section 7. 

Let 0,\. k ^ r, be the operating cost of the machine A, h in period k. With reference to 
our previous paper [51, we note that 

OH, - tt," - S, h , + M* 
and 

O** - - S? k + M, h k 
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where nj' is the purchase price (il is assumed thai the new machine .4/' is purchased ai the 
beginning of period /) and Sl\ is the salvage value at the end of period A of the machine .4," 
The term Ml\ denotes all the costs, which are specific to the technology /;. of operating the 
machine .4/’ in period A-excluding the cost of depreciation (or loss in salvage value from period 
A - I to A). 

To obtain a forward algorithm, we must solve the multiple alternative replacement model 
for any given finite horizon T. Suppose the shop uses n machines during the interval < 1,T> 

Let t 2 . i„- t . /„ be their salvage times and let h\, h 2 . h„ be their technologies. 

Note that r, > 1 for i = l, 2, ... , (n - 1) and t„ = T, since we assume that the shop goes out 
of business at the end of the given finite horizon T. Define t 0 = 0. Also note that 

<p.q> = {p, p + I.?! with nonnegative integers pand q with q > p. We can now state 

the finite horizon problem as 

(U Minimize £ £ 0*LY* 


h\, h 2 . h„. 

In the next section, we develop a forward dynamic programming algorithm to solve (1). 

3. FORWARD ALGORITHM 


We start with defining the following terms: 

Purchase Point. A period t is defined to be a purchase point (or P-pointl if a machine is 
purchased at the beginning of period t. 

Regeneration Point A period t is a regeneration point (or P-point) if a machine is sal¬ 
vaged at the end of period t. 

It is easy to see that any P-point is immediately preceded by an P-point. This property 
allows us to develop an efficient forward algorithm (2]. 

We need to introduce the following notation: 


cm 

C h (j, T) 


C ; "( T) 


(3) 

Cj(T) 


“ minimum cost of the T-period problem (1), 

= the total operating cost for the machine 
Af+i in the interval <j + l,f> 

= minimum total cost of the T-period problem when (/ + I) is the last 
P-point and is the technology of the last machine 

= C0) + C h (JT) 

= minimum total cost of the T-period problem with 0 + 1) as the last 
purchase period. 
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Clearly, 

(4) C,(7) = min C/(7) 
and 

(5) C(7) = min CAT) 

jt< o.T-i> ' 


Equations (2)-(5) complete the statement of the forward algorithm. It is possible to sim¬ 
plify the computation of C*(7); using (2) in (3) gives the following sequential procedure 


( 6 ) 


C,*< T) 


C*(r- 1) + 0; +1 .r for j < T- 1 

c(r- i) + o£ r y-r-r 


4. PLANNING HORIZON RESULTS 

To obtain planning horizon results, we need to assume that every technology is an 
improving technology. Thus, for technology A, we assume that 

(7) 0,-u > 0* k for k > j + \. 

In other words, we assume that the operating cost in a period for a machine of technology h is 
lower than the operating cost for an older machine of the same technology, except perhaps in 
the first period of operation of the new machine when the amount of depreciation on the new 
machine is usually high. 

Lower Bound Monotonicity Property: 

Let j*(T) denote the latest next to the last /?-point in the optimal solution of the 7- 
period subproblem; note that period T is the last /{-point. By the lower bound monotonicity pro¬ 
perty we mean that the lower bound of j*(T) increases monotonically with T. To prove this 
property, we must also define j* h (T) to be the latest next to the last /{-point for the optimal 
solution of the C h ( D-problem, where the C*! 7)-prob)em is the 7-period problem subject to 
the constraint that the last machine be machine of technology h. 

THEOREM I; Under the assumption (7) of improving technology over time, j* h (T) 
satisfies the regeneration monotonicity property, i.e., 

(8) y**(7 + 1) > j* h (T). 

PROOF; Let j* h (T + 1) - a and y**(7) - b. It is easy to see that C£(D ^ C£(D. 
Furthermore, if a < A then OjVi.r+i ^ 0j +1 . r+ |. 

If a < b, then using (6) we can write 

C*(7+ 1) - C*(7) + O a \,. r+l 
^ C»(7) + 0»\,. r+1 
- Q*(7+ 1). 
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This contradicts the optimality of the CjUT + l)-solution, implying that a > b. This completes 
the proof. 

THEOREM 2: If we define 
(9) X(D « min (/**(D}, 

then X(7~) represents a lower bound of j*(T + 1) and X(D is monotonically increasing in T. 
i.e., 

no) r(T + i) ^ x(D > x(r — i). 

PROOF: The proof easily follows from (8) and (9). 

Regeneration Set: 

The regeneration set of the C(D -problem can now be defined as 

(11) HIT) - set - |k(D, x(D + 1. T-l}. 

We note that a set is known as an R( D-set if it contains at least one R-point of an optimal 
solution to any problem with horizon T or longer. 

In many cases, it is possible to reduce the size of the /?(D-set. This reduction makes it 
easier to obtain planning horizons as can be seen in the planning horizon theorem below. The 
reduction procedure requires computations of the S'TD-sets for all h, where an S''( D-set is 
defined to be a set of regeneration points such that j* h (T + /) € S h ( D-set for all periods 
(T + /) with j* h i T + I) < T and / > 0. It is convenient to express this procedure by the flow 
chart on the next page. 

The reduced R (D-set can be obtained as follows: 

(12) R (D -set = U SH D-set. 

It should be noted that the /?(D-set in (12) cannot be bigger than the R( D-set in (11). We 
now prove the following important result for the set S found in the flow chart. 

THEOREM 3: If a period r $ the S-set found in the flow chart for technology h , then 
j*h{ t + i) tf or any value of / > 0. 

PROOF: If i ? the S-set found in the above flow chart, then there is a period u, 
T - 1 > u > t, such that C*(D < C,*(D. From the improving technology assumption, we 
also have 0,'' +u . < 0,+ i.* for k > T + 1. We can now write for any period T + /: 

T+l 

OUT + I ) - CUT) + £ 0* +l ,* 

T+l 

< C*(D + £ 0,* +u 

*-r+i 


- c,*(r + /), 

implying that j* h (T + /) & t. This completes the proof. 
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Planning Horizon Theorem: 

Let f*(r) denote the period when the first machine is salvaged in an optimal C(re¬ 
solution. We now state the planning horizon theorem which gives us a stopping rule for the 
forward algorithm to find the optimal replacement time in an infinite horizon problem. 

THEOREM 4 (Planning Horizon Theorem): If ,f*(r) = rf- constant for all r € R(T)- 
set (12) for some T such that j*(T) ^ 0, then the salvage time of the first machine in an 
optimal infinite horizon solution is t*. Furthermore, the optimal technology It* is determined 
by the condition 

03) cS f (/f> = c„or>. 


The first replacement time rf is called the planning horizon and T is called the forecast 
horizon. We need only to forecast the relevant technology up to period T to find the optimal 
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5. NUMERICAL EXAMPLE 


We now illustrate the planning horizon procedure by solving a numerical example This 
example assumes that machines of two different technologies are available for replacement at 
any time. The operating costs O,{ and 0,f are as shown below: 

O,) Or, 


2 _ _ J_ 


160 120 130 140 150 160 


300 60 70 80 90 100 

160 120 130 140 150 


300 60 70 80 90 

140 105 115 125 

/ 

200 50 60 70 

140 105 115 


400 30 40 

120 80 


400 30 

120 


2001 


Note that technology 1 can be interpreted as a labor intensive technology and technology 
2 can be interpreted as a capital intensive technology. The calculations of C,'(T ), C?(T), and 
C,(D are shown below: 



Sample calculations for T = 3 are shown below: 

C(j (3) - Of, + Ofj + Of 3 - 160 + 120 + 130 = 410, 
Cf (3) — C(l) + Of 2 + Of 3 = 160 + 160 + 120 = 440, 
Cf (3) = C(2) + Of 3 = 280 + 140 = 420, 
j*'( 3) - 0. 


NAVAL RESEARCH LOGISTICS QUARTERLY 


VOL 29, NO 3. SEPTEMBER 1982 






1*1 SWING 


SI M 11 


491 


C's (3) = Or, + Of : + O, 2 , = 300 + 60 + 70 = 430. 

C'f (3) = CU) + Of, + 0 2 2 , =• 160 4- 300 + 60 - 520, 

Cl (3) = C<2) + Of, « 280 + 200 = 480. 

/*'<3) = 0. A(3) = mini./* 1 (3). y* 2 (3)l = 0. 

C',,(31 = minfC'o (3), C’,?(3>) = 410. 

0,(31 = min(0,‘ (3). C, 2 (3)> = 440, 

0 ; (3) = min(0: (3). Cl (3)) = 420. 

./'*(3) = 0. /*<3) = 3. 

The planning horizon theorem is not satisfied for T < 5. For example, for T = 5. we 
have S 1 (53 = (4). S : ( 5) = (2,3.4), and /?(5) = {2.3,4). Since ,/*(2) * /*(3) * /*(5). the 
planning horizon theorem is not satisfied. For T — 6, we have s'(6) = (5), S 2 (6) — (2.5). and 

R (5) = {2,5). /*(2) = /*(5) = 2, so the planning horizon is 2 periods and the forecast hor¬ 

izon is 6 periods. 

6. RELAXING THE IMPROVING TECHNOLOGY ASSUMPTION 

Let /(/;) ^ 1 be a number associated with technology A such that 
0,'L|., + * 3? 0* J+k for k > t(h). 

This condition is likely to be satisfied by most improving technologies. The following 
result holds for this assumption: 

THEOREM 5: Let /*"( T) < T- tlh) (0 if T < /(/»)) be the latest period such that 
C, C,(T) for all j € <0, T - t(h)>, then j* k l T + 1) 3s /**( T). 

PROOF Assume to the contrary that i* h (T + 1) = a < i* h (T). It is easy to see that 
i > O," *h(T) + 1. T+ l and C^T) > C* Vt, r >m. We can now write 

C* *A, ri (r + I) - C* *h( Tt (T) + O ,* • /iin + |,r+i 
< C*m + 0* +l . r *, 

- C „(T + 1). 

implying that y**(T + 1) ^ /•*(D. 

Regeneration Set. 

Let A(D -* min [/**(7\>1, then the regeneration set of the C< D-problem can be defined 
as 

RID-set - |X(T), A(D + 1.T - 1). 

As in Section 4, it is possible to find a reduced regeneration set. For this we first find the 
S h <T )-set as below. The reduced regeneration set can be found by using (13). 
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Flow Chan lor Computing the S*( 74-set 


7. EXTENSIONS AND CONCLUDING REMARKS 

In this paper we have solved the infinite horizon replacement model with multiple alterna¬ 
tives. We have obtained planning and forecast horizons for this model under reasonably gen¬ 
eral conditions. Furthermore, the model can be easily adapted to situations where the set of 
alternative technologies is changing over time. This is done by setting O';, = =» for the ima¬ 
ginary A* machine for all./' < r for a technology which appears at lime / for the first time. If a 
technology h disappears in period i. then we let O' 1 , = °° for fictitious machines .4'’, for all 
./ > i. With these definitions, the model developed in the paper is applicable to the situation of 
changing sets of technologies. 

Another situation to which the model can be easily adapted is the situation in which a 
switch over cost k ah occurs whenever a machine of technology a is replaced with a machine of 
technology h. For this, we do require the assumption that k uh is separable, i.e.. 

k uh = k“ + k h . 

In this case, we can adjust the salvage value of the existing machine downward by amount k" 
and adjust the price of the new machine upward by amount k h . 

Finally, we must slate that the extensions of the machine replacement model m ..m r. . 
ous paper [5| can also be solved in the multiple .iltem.iiice case in a t.mU s;r.i.i:ti:t « 
manner 
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I. Adiri and D. Pohoryles 

TECHNION-lsrael Institute of Technology 
Hafa, Israel 

ABSTRACT 

This paper deals with flowshop/sum of completion times scheduling prob¬ 
lems, working under a "no-idle" or a "no-wait" constraint, the former prescribes 
for the machines to work continuously without idle intervals and the latter for 
the jobs to be processed continuously without waiting times between consecu¬ 
tive machines. Under either of the constraints the problem is unary NP- 
Complete for two machines. We prove some properties of the optimal 
schedule for n/2/F, no-idle/XC ( . For n/m/P, no-idle/XC, and n/m/P, no- 
wait/XC, with an increasing or decreasing series of dominating machines, we 
prove theorems that are the basis for polynomial bounded algorithms. All 
theorems are demonstrated numerically. 


INTRODUCTION 

The nonpreemptive, flowshop, sum of completion-times scheduling problem is: n jobs 

(Ji,J 2 . J„) have to be processed by m machines ... , M m ). Job J h 

/'“ 1,2. rt, consists of, at most, m operations (0,1,0 (2 .0 lm ). Operation 0,j which 

precedes 0 lj+l , has to be processed uninterrupted for t u time units, on Afj, j - 1,2. m. t 0 

is a nonnegative integer — t,j « 0 if 0 ,j is missing and positive if it exists.* Two operations of 
the same job cannot be processed simultaneously and a machine may process at most one job at 
a time. Find the operation sequence on each machine, that obeys the problem constraints and 
minimizes the sum of completeion times. The problem is designated n/m/F/ZC h t where F 
stands for flowshop discipline and C, is the completion time of job 

The no-idle constraint: machine M k , k — 1,2, .... m, works continuously without idle 
intervals. 

The no-wait constraint: job /,,/■■ 1,2, ... , n is processed continuously without waiting times 
between consecutive machines. 

Both constraints arise in real life situations. Examples of such scheduling are: (i) under a 
no-idle constraint — use of very expensive equipment (.e.g., a computer and its peripheral dev¬ 
ices) with the fee determined by the actual time consumption; (ii) under a no-wait constraint 
— in metal-processing industries (e.g., hot rolling) where delays between operations interfere 
with the technological process (e.g., cooling in the above case). 


•For further elaboration of the meanings and influence of zero processing times see Hefetz and Adiri (41. 
tFor notation and classification of scheduling problems we follow Lenstra ($1 and Rinnooy-Kan (6). 
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While the no-wait case has been widely discussed in the professional literature, [11 [61, 
the no-idle case is defined here for the first time. This paper deals with flowshop disciplines 
with a no-idle or a no-wait constraint where the objective function being minimum sum of 
completion times. 

Garey, Johnson and Sethi [31 proved that n/2/F/ZC, is unary NP-Complete by construct¬ 
ing an instance that is no-idle and no-wait, and thus proved at the same time that of n/2/F, 
no-idle/XC, and n/2/F, no-wait/XC,. Section 2 deals with some properties of the optimal 
schedule for n/2/F, no-idle/XC,. Sections 3 and 4 are devoted to proofs of theorems that 
underlie polynomial bounded algorithms for n/m/P , no-idle/XC, and n/m/P, no-wait/XC, with 
an increasing or decreasing series of dominating machines, respectively. 

1. COMPLEXITY OF FLOWSHOP/NO-IDLE OR NO-WAIT/SUM OF 

COMPLETION TIMES PROBLEMS 

Garey, Johnson and Sethi [3] proved that 3-partition is reducible to n/2/F/Z.C,. How¬ 
ever, as the constructed instance of n/2/F/XC, happens to be no-idle and no-wait, thus at the 
same time proved the unary NP-Completeness of n/2/F/ZC it n/2/F, no-idle/XC, and n/2/F, 
no-wait/XC,. Moreover, only minor modification of the proof is needed for proving the NP- 
Completeness of the first two problems where missing operations are prohibited. Specifically, 
replacement of zero processing times on A/, (missing operations) by c > 0 (existing operations 
with infinitesimal processing times) and shifting of all « to the beginning of the schedule on 
A/], proves the unary NP-Completeness of n/2/F, t u > 0/XC, and n/2/F, t u > 0, no-idle/XC,. 
The complexity of n/m/F, t,j > 0, no-wait/XC, for fixed m > 2 is an open problem. 

2. CONSTRUCTION OF A NO-IDLE SCHEDULE 

We distinguish two ways of constructing a no-idle schedule. Both are implemented con¬ 
secutively on the machines, starting with the second, proceeding to the third and so on until 
the last machine. 

(0 Right shifts. We shift to the right every operation that precedes an idle interval 
(maintaining the constraints of the problem on M k , k > 2, further right shifts might be 
needed), until all idle intervals have been eliminated. 

(ii) Left shifts. At first we schedule the operations on M k , k > 2 in a single stretch (a 
block) without idle intervals, starting when the last job on A/*_| has been completed. After¬ 
wards we apply a maximum left shift (without violating the constraints) to the whole block. 

Let us define a blocking job on M k , k ^ 2, in a no-idle schedule as the first job on M k 
that prohibits shifting to the left. 

3. PROPERTIIES OF THE OPTIMAL SCHEDULE FOR nf2jF, NO-IDLE/XC, 

Conway, Maxwell and Miller [21 fundamental theorem for flowshop scheduling that states 
that an optimal schedule exists for a n/m/F/b problem (8—any regular measure of perfor¬ 
mance) with the same processing order on the first two machines [2, p. 81], holds for the case 
under discussion. Thus, the set of permutation schedules for n/2/F, no-idle/XC, is a dominant 
one. 
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THEOREM 1: If in the optimal schedule for n/2/F, no-idle/EC,, the blocking job, is the 
last one, then the optimal schedule is according to SPT (Shortest Processing Time) on M 2 , 
except for the blocking job that is the one with the minimum processing time on Af 2 . 


PROOF: Let us denote, a, - b, - r (2 ; S lk starting time of J, on A/*; and the square 
brackets indicating the place of a job in the sequence, for example, S[n 2 —starting time of the 
first job on A/ 2 . Since the blocking job is the last one, we have (see Figure 1). 

M %2 " £ «(/! - L *l<! ” K + *l«l- 

where 

K — jT a, - b, — Constant. 




w i 


S |112 


3 In] 2 


The sum of completion times takes the form 

(2) C, - S |I12 + b u] + S (112 + + ft l2t + ... + S, 112 + b[,\ 

” "$1112 + £(»-/ + 1 ) 6|/|. 

Substitution of (1) yields 

(3) C) “ nK + nA[ B | + (n — i + l)A(,| 

■ nK + (« + 1)A|„| + (n — / 4- l)/>[ ( ). 

To minimize (3) b\„\ should be the smallest and 6(<], / — 1.2. n — 1, a nondecreasing 

sequence, thus the optimal schedule is given by the sequence 

ft l»l < *lil < *121 < ■ • • < *|»—il- LI 


THEOREM 2: The jobs that (i) precede (ii) succeed the blocking job in the optimal 
schedule for n/2/F, no-idle/XC, are ordered according to SPT on Af 2 , provided the no-idle con¬ 
straint is not violated. 

PROOF: Let us assume that the blocking job in the optimal schedule is the w-th in the 
sequence. 
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M, Jmi 

m» 1 

Figure 2. Optimal schedule for n/2/F, no-idle/IC, where the blocking job is J l# |. 


(i> -^1112 “ X a U\ ~ 5J *1*1- 

The situation here is rather similar to that of Theorem l with w replacing n. Thus, in the 
optimal schedule we have 6m < 6m < ... < 6 ()r _ 1 j. 

(ii) For / - ir + 1, n + 2, , n, we have 

C (l]"C(ir|+ X *1/1* 
y-»+l 

Thus, 

X C W) “ ~ w ) C lirl + X ~ ’ + 1 ^*I*1* 

/-ir+I i-ir+1 

We have that in the optimal schedule the jobs that succeed the blocking job obey the SPT rule 
on A/ 2 , provided the no-idle constraint is not violated, 6 lw+II < 6[ jr+2 ] < ... < b\ n \. □ 

4. n/m/P, NO-IDLE/IC,, m > 2, WITH A SERIES OF DOMINATING MACHINES 

Machine M k dominates M, if 
(4) min t ik > max 

In abbreviated notation 

M k > M r . 

A scheduling problem with a series of increasing [decreasing] dominating machines is one 
where M, <M 2 <••• <M m lM\ > M 2 > ...> A/J. 


For n/m/F, no-idle/IC,, m > 2, the set of all permutation schedules is not a dominate 
set. However, for the sake of solvability (development of polynomial bounded algorithms for 
special cases) we confine our search for optimality to permutation schedules, and the problem is 
designated n/m/P, no-idle/IC,. 

Note that for m - 2 the set of permutation schedules is a dominate set. 

Let y lk be the processing time of the Mh job on M k where the order is according to SPT 
on the last machine, M m . 
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THEOREM 3: The optimal permutation schedule for n/m/P, no-idle/IC, with an 
increasing series of dominating machines* is according to SPT on the last machine, except for 
the first job that satisfies 

Min | Vj - n £ y jk + (j - l)y, m - y, m J. 

PROOF: For this case we have 
(5) J C, - nS [Um + (« - / + 

where 5[ 1]m is the starting time of the first job on the last machine, and 
M S[ilm “ 'ill/- 

Thus, for a given first job the optimal permutation schedule is SPT on M m , 

'l2)m < '[31m < • • • < 'Mm- 

Selecting the >th job in the SPT sequence on M m to be the first, we have 

(?) 1 - » I yj> + O' - 1 )v Jm - I Y*» + |(—i+ l)y,m. 

Since the last term in (7) is not affected by the choice of the first job, the latter is taken so as 
to satisfy 

m ! n | Vj *» j Vjk + 0- Dy>* “ £ T/mj O 

EXAMPLE 1: A 5/4/P, no-idle/IC, with an increasing series of dominating machines 
with processing times as per Table 1. 


TABLE i — Processing Times 


\^Jobs 

Machines\^ 

J\ 

h 

h 

J 4 

Js 

M\ 

3 

1 

3 

2 

2 

Mi 

4 

5 

6 

4 

6 

Mi 

7 

9 

8 

8 

9 

__1 

11 

13 

14 

10 

12 


M 4 > Mi > M 2 > M u thus we have an increasing series of dominating machines. 
SPT sequence on M 4 is 4-1-5-2-3. 


’Since M k+I > M k , k — 1,2, .... m - 1, the no-idle constraint is not effective and for the case under discussion, the 
two problems n/m/P, no-idle/IC, and n/m/P/ZC, are equivalent. 
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The first job, say J n satisfies V, — Min V r is the Ath in the SPT sequence on A/ 4 . 


V\ - 5(2+44-8) . - 70 

V 2 - 5(3+4 + 7) + 11-10 - 71 

Vy - 5(2+6+9) + 212- (11 + 12) -88 

— 50+5+9) + 3*13— (10+11 + 12) -81 


V s - 5(3 + 8+8) + 4 14 — (10+11 + 12 + 13) - 95 


The job chosen to be the first is the first in SPT sequence on (F, - min V } ), namely J k . 
The other jobs are scheduled according to SPT on A/ 4 . Thus, the optimal permutations 
schedule is 4-1-5-2-3 with T C, - 240, and takes the form as per Figure 3. 


*rns0ii 


UM 1 5 | 2 | 3 | 

<*■> 

CM 

in 


n | ’ | s i 

1 2 ] 2 | 

24 35 4 

Figure 3. Optimal schedule for Example 1 

7 60 74 


THEOREM 4: The optimal permutation schedule for n/m/P ’, no-idle/IC, with a decreas¬ 
ing series of dominating machines is according to SPT on the last machine except for the last 
job that satisifies 

nunj V} - « y Jk -in- j)y Jm + y, m |. 

PROOF: Since M\ > M 2 > ... > M m the blocking job on M k , k - 2,3.m, is the last 

one, we have 

(8) S[il m - Jj ($li)<*+i) ~ Sill*) 

" Z |l 'i/i* ~ X 'i/n*+i)J 


-I 


(T k - T k 


+ '(»l(*+t)) • K + 2} 'l»l(*+i)’ 


where T k is the total processing time demanded on M k , T k - constant - and 

X- 7-,- T n . 
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Substitution of (8) in (S) yields 

X Cj - nK + n £ /(„K* + i) + T (« - / + 

i-l *«l 

Given the job in the last (n-th) place, we have 

go) Jj c ( - a:, + j (.-/ + i)/ f(lw , 

where Afi = constant - nK + n £ /foiu+n + /(„)„. 

Minimum of (10) is obtained by an SPT sequence on M m , 

'lllm < 't21m < ••• < 'lit-llm- 

Choosing the >th job in the SPT sequence on M m for the last place, we have 
D X C, - nK + n T y j(k + „ -in- j)y Jm + £ y im + £ (n - i + l)-y, m . 

i~ *-l i-j +1 ” 

The last term on the right hand side of (11) is independent of the choice for the last place, thus 
the job is taken so as to minimize 

“ " X ?/<*+» “ “ Jbjm + X, D 

COROLLARY 1: The optimal schedule for n/2/F, no-idle/IC, where M\ > M 2 is 
according to SPT on M 2 except for the last job that is the one with minimum processing time 
on Mi- 


PROOF: As was pointed out the two problems n/2/P , no-idle/IC, and n/2/F, no- 
idle/IC,, are equivalent. Substitution of m « 2 in Theorem 4 yields 

v \ ” "I'" I v i “ Jyj2 + X 

' [ i-J +1 

Thus, the last job is the first in SPT sequence on M 2 . This result is in agreement with 
Theorem 1—since Af, > M 2 the blocking job is the last one and should be the smallest. □ 

EXAMPLE 2: A 5/4/P, no-idle/IC, with a decreasing series of dominating machines 
with processing times as per Table 2. 

TABLE 2 — Processing Times 


y, h h j 4 j 5 

11 13 14 10 12 

7 9 8 8 9 

4 5 6 4 6 

3 13 2 2 
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Mi > M 2 > Mj > M 4 , thus we have a decreasing series of dominating machines. 
SPT sequence on A/ 4 is 2-4-5-1-3. 

The last job, say J h satisfies V t — Min Vj, is the /-th in the SPT sequence on M A . 
V, = 5(9+5 + l) - 41 + 2+2+3+3 = 81 


V 2 = 5(8+44-2) - 3-2 + 2+3+3 = 72 

F 3 = 5(9+6+2) - 2-2 + 3 + 3 - 83 

y 4 = 5(7+4+3) —1-3+3 =70 

V s = 5(8+6+3) = 85 


The job chosen to be the last is the fourth in SPT sequence on Af 4 , ( V 4 = min Vj), namely J\ 
The other jobs are scheduled according to SPT on A/ 4 . 


Thus, the optimal permutation schedule is 2-4-5-3-1 with C, = 343 and takes the form a< 
per Figure 4. 


I a 1 s | 


X=L 


m. 


i 5 i 3 m 


2 U 1 S I 3 I 1 j 


_§ 


imn 


646668711 


Figure 4. Optimal schedule for Example 2. 


5. n/m/P, NO-WAIT/IC,, m > 2, WITH A SERIES OF DOMINATING MACHINES 

We recall that "no-idle" constraint prescribes for the machines to work continuously 
without idle intervals, while the no-wait constraint prescribes for the jobs to be processed con¬ 
tinuously without waiting times between consecutive machines. 

The set of all permutation schedules for an "F, no-wait" problem (fy > 0, missing operations 
are allowed) is not a dominant set. However, we confine our search for optimality to permuta¬ 
tion schedules, namely, the problems under discussion in this section, are particular cases of 
n/m/P, no-wait/EC,. Note that for t tJ > 0 (missing operations are not allowed) a feasible solu¬ 
tion for an " F no-wait" problem is a permutation one, thus "F, t,j > 0, no-wait" and "P, t 0 > 0, 
no-wait" are equivalent. As was previously pointed out, n/2/F, no-wait/EC, is NP-Complete 
while the compexity of n/m/F, no-wait, t y > 0/EC,, m > 2, is an open problem. 
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THEOREM 5: For an increasing series of dominating machines the two problems n/m/P, 
no-wait/IC, and n/m/P, no-idle/IC, have the same optimal sequence. 

PROOF: For n/m/P, no-wait/IC, with an increasing series of dominating machines, we 

have 

C[,| * £ f(in + 21 't/Iio- 

It b readily shown that (12) leads to (5) and (6). □ 

A direct consequence of Theorems 5 and 3 is that the optimal schedule for n/m/P, no- 
wait/IC, with an increasing series of dominating machines is SPT on M m except for the first 
job that satisfies 

Min j Vj - n £ Ty* + </ - »>ym ~ Jj y, m |- 


EXAMPLE 3: A 5/4/P, no-wait/IC, with an increasing series of dominating machines 
with processing times as per Table 1 (Example 1). 


* 10 . 


HUE 


a 


HL 


mm 


m 


m 


*3_ 1 * I I \ H 5 I I l 1 [ 3 | 

•«l-1 4 I 1 I I I 2 1 

24 35 47 60 

Figure 5. Optimal schedule for Example 3. 


V_J 


The calculations and the optimal sequence are the same as in Example 1, but the resulting 
schedule is as per Figure 5 (compare with Figure 3). 

THEOREM 6: The optimal schedule for n/m/P, no-wait/IC,, m > 2, with a decreasing 
series of dominating machines* is according to SPT on M x . 


•Since M k > M t+I , * - 1,2, .... m - 1, the no-wait constraint is not effective ; for the case •• ir discussion the 
two problems n/m/P, no-wait/IC, and n/m/P/XC, are equivalent. 
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PROOF: Since M k > M k +\, k « 1,2. m - l, we have for the optimal schedule 

(13) £ C, - filly + '|1J1 + £ '12I> + • • • + £ Hill + £ '(»ly 

“ 5^ 2^ 'I'U + £ ~ 

The first term on the right-hand side of (13) is a constant, and the sequence 
'(ill ^ 'mi < < ' 1 ,,-ui minimizes the second term. Morever, r^n > otherwise if 

/(„li ^ max f(,n the value of the second term can be reduced by interchanging the last job with 

4h- 'L/li “ m f x 'l/ii- a 

EXAMPLE 4: A 5/4/P, no-wait/EC, with a decreasing series of dominating machines 
with processing times as per Table 2. Since M t > M 2 > M 2 > A/ 4 we have a decreasing series 
of dominating machines and the optimal sequence is according to SPT on M t , 4-1 -5-2-3, with 
jr C, ■» 247. The optimal schedule takes the form as per Figure 6. 
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ABSTRACT 

For solving transportation problems essentially three types of methods are 
known: primal methods, the Hungarian method and the shortest augmenting 
path method. In this paper we present the specialization of these approaches to 
the bottleneck transportation problem and report some computational experi- 


1. INTRODUCTION 

The classical transportation problem (TP) can be formulated in the following way. There 
are m supply points and n demand points with supply point / capable of supplying amount 
<?,(/ - 1. m) and demand point j having demand bj(J - 1. n) with £ a, - £ b } . 

Find the least cost transportation pattern from the supply points to the demand points when the 
cost of transporting a unit amount from / to j is c u% i.e.. 



min L L C ij ’ X u 

subject to 


(1.1) 


(/' — 1. m ) 


(1.2) 


O' ” 1. n) 


(1.3) 

x,j > 0 

(»- 1. m\j- 1. . 

.. . n] 


A related problem is the bottleneck transportation problem (BTP). Here a transportation 
time t,j is specified between each supply point /' and each demand point j. Now it is required to 
find a transportation pattern which minimizes the total time necessary for transporting the 
goods from the supply to the demand points. Thus, 

(1.4) min max {t u \x tJ > 0} subject to (1.1), (1.2) and (1.3) 

Problem (1.4) is often called "time-transportation problem." It occurs in connection with tran¬ 
sportation of perishable goods, with the delivery of emergency supplies or when military units 
are to be sent from their bases to the front. 
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For solving TP we know essentially three types of methods 

— primal methods 

— the "Hungarian" method 

— the shortest augmenting path method. 

Nearly every approach for solving TP can be modified to treat BTP. 

So for instance the primal method of Hammer [10], [lOal, for solving BTP is strongly 
related to the method of Klein [11] for solving TP. Starting from a feasible solution both 
methods are looking for "negative cycles" to improve the actual solution. 

In this paper we present these three basic approaches for solving BTP and we report some 
computational experience with these methods. 

2. THE HUNGARIAN METHOD 

This method for solving BTP was first proposed by Garfinkel and Rao [8], Inspired by a 
work of Edmonds and Fulkerson [51 on general bottleneck problems, they called the procedure 
the Threshold Method. The algorithm can be described in the following way: First determine a 
"good” lower bound z for the optimal objective value z. Then define a network Jf{z) in the fol¬ 
lowing way: Let V - {s,/}0 M 0 N with A/ — {1,2, .... m } and N = {m + 1. m + n). 

be the nodeset of Jf(z). Now every node i € M is connected with s by an arc, respectively, 
each m + j € N with /. A pair (i,m + J) is connected by an arc only if r y < z holds. 

The arc-capacities d ki i are defined by 

d yi :- a* (/' - 1.m) 

dt. m+J 00 d - 1. m\j- 1. n) 

d m+J , bj 0 -1.«)• 

Now a maximum flow / - (f 0 ) from source s to sink ris determined using the labeling method 
of Ford and Fulkerson [6], If for the maximal flow value v - £ a, holds, the optimal solution 

for the BTP is obtained defining 

x u :— (/' — 1. m\j — 1.n) 

If v < I a, the labeling procedure yields simultaneously a minimal cut (X.X) in Jfiz) 
with s € X and / 6 X. This cut has the property 

f y 2 -► (i,m + j) € (X,X) 

To obtain a solution to (l.l)-(l.3) it is therefore necessary to use an arc (t',m + j) with t 0 > z 
Hence, we determine 

z* :— min {t,j\i € X, m + j €X) > z 
Define i :- z* and repeat the process. 

After (m ■ n) iterations, at most, an optimal solution is obtained. The algorithm is sum¬ 
marized in the following flow-chart: 
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Figure 1. Flow-chart of the Hungarian method 


During the algorithm the network Jf(z) need not be constructed explicitly. All operations of 
the labeling procedure can be performed directly on the "time" — matrix T - U u ). 

3. A PRIMAL METHOD 

The first method for solving BTP proposed by Barsov [1] is a primal algorithm which 
proceeds "dual" to the Hungarian method in some sense. In the course of the Hungarian 
method the "threshold" z is successively increased until a feasible solution can be found. In this 
method a feasible solution is always at hand and the threshold-value is successively decreased 
until no further solution can be found. Starting from a feasible solution x - (x fJ ) with objec¬ 
tive value z a cost-matrix D » (dy) is defined by 

| 0 if ty < z 
d » | 1 else 

Now a TP with cost-matrix D is solved. If the optimal solution y «■* (y 0 ) has an optimal 
value z(y) > 0 the solution x — (x y ) is optimal for BTP. Otherwise 

z* max {r y |> y > 0} < z holds. 

Define z z* and x y and repeat the process. 

Computational experience shows (8) that this method is inferior to the Hungarian 
method. Recently, Finke and Smith (7) developed an improved primal method. 
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Starting from a lower bound z and a "good" start solution x - (x u ), a cost-matrix 
D - (d u ) is defined by 

djj :« | — ^ y j with W greatest integer < x. 

Now the associated TP is solved using the well known MODI-method and x - (xy) as 
start solution. If the optimal solution y — (y 0 ) has objective value z(y) « 0, then y is optimal 

for BTP. Otherwise, we consider the dual variables / — I, .... m and v y , j — 1. n 

associated with the optimal solution y. Then 

u, + \j < dij i = 1, ... , nr, j — 1. n 

y,j > o -*• u, + \j - d,j 

and 

z ’min {/<, !«, + v y > 0) > z holds. 

It can be shown that z' is a lower bound for the optimal value for the BTP. Define z :« z' 
and repeat the process. In an implementation the actual basic solution of the TP is stored as a 
tree using the list structures proposed by Glover and Klingman [9], Srinivasan and Thompson 
[12]. Due to this technique primal methods are superior to other methods in the case of TP. 

The following flow-chart summarizes the algorithm of Finke and Smith. 



Figure 2. Flow-chart or the primal method 
of Finke and Smith 
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4. THE SHORTEST AUGMENTING PATH METHOD 

This method was proposed by Totnizawa 114] for solving TP and applied to BTP by Derigs 
and Zimmermann [4], A theoretical foundation of this method can be found in [2], During 
the procedure so-called /c-subproblems (P k ) have to be solved for 1 ^ k < nr. 

(P k ) min max {fy \x u > 0| 

I *'" 1 ’ 2 .* 

fa x o \ 0 else 

Jjxy < bj j - 1. n 

Xjj'Z 0 i - 1.m ; j - 1. n. 

The solution for (P t ) is obvious. Starting from an optimal solution x - (x„) for (P k ) 
with k — 1, , m - 1 an optimal solution for (P* + i> is obtained by means of augmenting 

paths. For this purpose we partition (1,2, .... n] *» /V, 0 N 2 with 

N t = |y I Xjj - ftyj "saturated columns.” 

Now a sequence Ax of mutually distinct matrix-entries 

(k + OWi), (/Wj), U 2 J 2 ) .(»,J r+ i) 

is called an augmenting path with respect to x - (x j; ) if 

x^ j" > 0 for q - 1.2. r and 

j q € N, for q «■ 1,2. r and y.+i € N 2 . 

The capacity of Ax is defined by 

cap (Ax) min — JJxy r+1> e) with 

e min (x/ fj \q - 1,2.r}. 

Let us define x x © Ax by 

x^ — cap (Ax) for (ij) - (/| Ji).(4 Jr) 

Xjj + cap (Ax) for (/J) € Ax\l(i,J,)l 9 - 1, .... r) 

Xjj else 

and 

w (Ax) max {t u \x,j > 0). 

Let D k be the set of all augmenting paths with respect to x - (x tf ) and Ax € D k with 

(4.1) w(Ax) < w(Ax) for Ax € D k 

(4.2) cap (Ax) - a* +I 

then x x © Ax is an optimal solution for (P k+ ,). 
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An augmenting path Ax € D k with property (4.1) can be determined using a modified 
Dijkstra method. If Ax does not fulfill (4.2) we split (P k +\) into two subproblems (P k+ \) and 
(/*+ 1 ) with 0 * + , cap (Ax) and a k+] a k+l - a k+ ,. Thus, (P k +\) can be solved in a finite 
number of steps. 


Computational experience shows that the method can be improved with the aid of the fol¬ 
lowing starting procedure. First we determine a "good" lower bound z for the optimal value z 
and a "good" partial solution x - (x„), i.e. 


(4.3) Jj x tJ < a, 

(4.4) Jj x tj < bj 

(4.5) x tj > 0 
with the property 


, / — 1 , .... m 


, / - 1, ... , m, j - 1. n 


(4.6) x tJ > 0 -► t u < 2 

such that 

5S x u is large. 


A partial solution x with (4.3) • (4.6) is called ^-feasible. Computational tests have shown 
that it is not useful here to determine a z-feasible solution with maximal possible sum. Some 
simple heuristics will give solutions which are within 5% from optimality much faster and, with 
respect to the running time of the complete procedure, this is of advantage. The complete pro¬ 
cedure is summarized in the following flow chart. 
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5. LOWER BOUNDS AND STARTING PROCEDURES 

All three methods presented above start with a lower bound z for the optimal value z. 
Computational experience shows that the bound proposed by Garfinkel and Rao is favorable. 
Row-thesholds z, and column-thresholds z j are computed with 

z, min max [t,j\x lj > 0} subject to 



0 < Xy < bj for /' = I, ... , m 
and z j defined analogously. 

Then z max {z|,Z 2 , .... z„, z'.z 2 , .... z"} is a lower bound. This bound is easy to 
calculate and yields a good estimation for the optimal value in most cases. 


The primal method of Finke and Smith heeds a good starting solution x — (x y ). They 
propose the following method. First a z-feasible partial solution x ■ (x 0 ) is constructed. Such 
a partial solution is used for the shortest augmenting path method, too. A z-feasible solution 
can be obtained in the following way. For every row i (column j) the number T ( ( P) of feasi¬ 
ble matrix entries iij) subject to t 0 < z is computed. Now the row i 0 (column j 0 ) with 
minimal T l(j > 0 (T 70 > 0) and positive a i(j (b J() ) is determined. If no such row (column) exists 
the partial solution is completed. Otherwise, determine in row i 0 (column Jo ) the feasible entry 
('Wo) w >th minimal number P 0 (T, 0 ) and positive b Jo ( a i(> ). Then define 


\.Jo •- min {a, 0 . b Jo ) 


Update Tj and P and repeat the process. 

This way a good z-feasible partial solution Jf - (x u ) is determined. If x - (x iy ) fulfills 
(1.1) and (1.2) then it is an optimal solution for BTP. If the partial solution is not optimal in 
the second step of the starting procedure proposed by Finke and Smith, the remaining supply is 
distributed using a north-west-corner rule. Finke and Smith call their procedure threshold- 
totals-method. 


The following flow-chart summarizes the starting procedures for the three methods. 
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FtcuRE*4. Flow-chart of the surting procedures. 


6. COMPUTATIONAL RESULTS 

FORTRAN IV implementations of all three methods were tested on a CDC CYBER 76 of 
the Computer Center of the University of Cologne. We used the following codes: 

HUNGAR - a version of the Hungarian method which Professor Garfinkel kindly made 
available to us 

PRIMAL - a version of the primal method which Professor Finke kindly made available 
to us 

SAP - an improved version of the BTP-code which is listed in (4]. The improve¬ 
ment concerns mainly the starting procedure. 


The HUNGAR-code uses the maximal cardinality NL of expected postive x u 's in any 
column as input data. The running time of the program is highly dependent on the choice 
NL since the value and position of the x,/s are stored in NL x N matrices. Thus, the choice 
smaller "NL" decreases the running time. But if the actual number exceeds NL, termination 
occurs. Choosing NL ” M/2 all problems were solved. 


We considered rectangular as well as quadratic examples. The integer coefficients of the 
time matrix T, the supply vector a and the demand vector b were generated by a machine 
independent uniformly distributed pseudo random number generator in the interval ll,2 J1 -1]. 
Afterwards these numbers were transformed into the intervals [1,6] with 
b - 10, 100, 1000, 10000, 2 JI -1. Twenty-five examples were generated for each combina¬ 
tion to calculate the mean running time. Numerical experience showed that the running time 
of HUNGAR and SAP for solving rectangular problems with m > n is less than for the 
equivalent transposed problem with m < n. Therefore, we recommend that rectangular prob¬ 
lems always be solved in the form with m > n. 


NAVAL RESEARCH LOGISTICS QUARTERLY 


VOL. 29, NO. 3, SEPTEMBER 1912 


■ o o 












SOLVING BOTTLENECK TRANSPORTATION PROBLEMS 


513 


The following table shows the mean running time for quadratic problems with 
m - n - 100. 


TABLE 1 — Mean Running Time in CPU-Secondsfor (100 x 100) BTP 


a„bj 


1-10 

1-100 

1-1000 

1-10000 

l-<2 31 -1) 


SAP 

.112 

.177 

.177 

.192 

.181 

— 

PRIMAL 

.106 

.124 

.124 

.133 

.126 


HUNGAR 

.226 

.400 

.383 

.478 

.504 

g 

SAP 

.107 

.120 

.144 

.171 

.150 


PRIMAL 

.094 

.106 

.105 

.106 

111 


HUNGAR 

.332 

.582 

.611 

.588 

.672 

"s 

SAP 

.089 

.140 

.188 

.121 

.223 

o 

PRIMAL 

.092 

.098 

.119 

.106 

.114 

- 

HUNGAR 

_333J 

.439 

.548 

.603 

.713 


The mean running time of PRIMAL is significantly better than the other's. For SAP this 
is simply caused by a single "ill conditioned" problem for which the running time is more than 
tenfold the running time of the other 24 examples of the group whose running time is compar¬ 
able to those of PRIMAL. 

SAP and PRIMAL are shown to be relatively insensitive to the range of the parameters 
di.bj and while the running time of HUNGAR is more than doubled when the range for r„ is 
increased from b - 10 to b - 2 31 -1. 

Then we modified the problems subject to 

t\j ■ f,t • 0 for / ■» 1. m and j — 1. n and 

a l - fe| - max (max {o ( J, max (ft,)). 

For these perturbated problems the starting procedure has no effect and thus the computational 
behavior of the "pure" algorithms is shown. 

Table 2 shows that PRIMAL and HUNGAR are more dependent on the quality of the 
starting procedure and the range of the time values r,, than SAP. 


TABLE 2 — Mean Running Time in CPU-Seconds for the Perturbated 
(40 x 40) Problems_ 


a,,bj 


1-10 

1-100 

1-1000 

1-10000 

l-(2 31 —1) 

0 

SAP 

.014 

.105 

.122 

.117 

.121 

~7 

PRIMAL 

.015 

.118 

.369 

.446 

.448 


HUNGAR 

.016 

.269 

.608 

.720 

.686 


SAP 

.014 

.153 

.196 

.202 

.169 


PRIMAL 

.014 

.133 

.389 

.483 

.504 

4 

HUNGAR 

.026 

.436 

.790 

.922 

.988 

9 

SAP 

.014 

.161 

.179 

.217 

.164 

© 

PRIMAL 

.015 

.132 

.398 

.483 

.506 


HUNGAR 

.025 

.422 

.834 

.937 

.983 
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Then we applied the algorithms to bottleneck assignment problems i.e., problems with 
a, - b, - 1. 

TABLE 3 — Mean Running Time in CPU-Seconds for (WO x 100) 

Bottleneck Assignment Problems 



, u 

1-10 

1-100 

1-1000 

1-10000 

1 - (2 31 -1) 

S E 

SAP 

.074 

.074 

.074 

.073 

.071 

E 5 

C £> 

00 O 

PRIMAL 

.083 

.092 

.100 

.096 

.105 

8 £ 

< 

HUNGAR 

.055 

.103 

.143 j 

.140 

.155 


For assignment problems SAP needs at most n/2 iterations while PRIMAL will perform 
a number of degenerate pivot operations. Again PRIMAL and SAP show a more robust nature 
with respect to the range of the time parameters ty than HUNGAR. 

More computational results can be found in [3] and (4]. 


7. CONCLUSIVE REMARKS 

In this paper we presented three basic methods for solving the bottleneck transportation 
problem. It is beyond the scope of this study to give an entire review of all the different 
methods which are available for solving BTP. Nevertheless, we think that the presented pro¬ 
cedures are representative in the sense that the main approaches for tackling combinatorial 
optimization problems are specialized to the bottleneck transportation problem. 

The computational tests indicate that PRIMAL and SAP outperform HUNGAR in gen- 
\ eral. The optimal choice between PRIMAL and SAP is data dependent and should be specially 

| made for every problem to be solved. If it can be expected that the start-heuristic is perform¬ 

ing well, PRIMAL is recommended. Otherwise, SAP is of advantage. 

We want to close our discussion by referring to another computational study. Werner 
, , [15] compared versions of SAP and PRIMAL on special problems the data of which was 

defined by some structured transportation networks incorporating time tables for transportation 
vehicles, waiting-times and different regions with different magnitudes of demand and supply. 
Werner reports that the running time for PRIMAL was highly dependent on the distribution of 
data and the more the data was not uniform the more SAP became superior. 
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SEQUENCING n JOBS ON TWO MACHINES WITH SETUP, 
PROCESSING AND REMOVAL TIMES SEPARATED 


Dileep R. Sule 
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ABSTRACT 

The paper extends the machine flow-shop scheduling problem by separating 
processing time into setup, processing and removal times. 


In a recent paper Yoshida and Hitomi [2] extended the classic two machine scheduling 
problem first introduced by Johnson (1], They allowed setup times to be independent of the 
processing times. This paper is the further extension of the Yoshida and Hitomi model; it 
allows for separation of processing time into setup time, processing time and removal time for 
each job on each machine. 

For example, consider a machine shop environment. Operations associated with each job 
on each machine when the machine/operator is available could be summarized as follows: 

1. Setup time independent of the unit to be processed. This operation consists of activi¬ 
ties such as obtaining the blueprints, procuring the necessary tools, fetching the 
required jigs and fixtures and setting them on the machine. 

2. Setup time that is a unit dependent. This operation includes the time required to set 
the unit in the jigs and fixtures and to adjust the tools as required. 

3. Processing time. 

4. Removal time dependent on the unit. This operation includes the times for activities 
such as disengaging the tools from the unit, and releasing the unit from the jigs and 
fixtures. 

5. Removal time independent of the unit. This operation includes activities such as dis¬ 
mantling the jigs, the fixtures and/or tools, inspecting/sharpening of the tools, 
returning them to the central depository, cleaning the machine and the adjacent area. 

Since activities 2, 3 and 4 are unit dependent, their times could be combined and desig¬ 
nated as the processing time. However, times for the activities 1 and 5 are independent of the 
unit to be processed. 


/ 
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Formulation of the Problem 

Let, 

P" - processing time of job /' on machine m. This includes times for the activities 2, 3 and 
4 described before, i — 1,2, ... #i, m — 1,2. 

S" — setup time independent of the unit, i.e., activity 1, for job / on machine m. 

RI" -removal time, i.e., activity 5, for job i on machine m. 

C” -completion time of job i on machine m. 

In addition, in the mathematical development, it is assumed that the jobs are designated 
so that the job in the Ah position of the processing sequence is job /. Since the setup is 
independent of the unit, if there exists an idle time on machine II, setup on machine II can be 
done before the unit is available from machine I. The unit is available from machine I when 
processing on the unit is completed; however, the machine is not available for the next job 
until the removal operation on the machine is finished. Figure 1 shows these activities graphi¬ 
cally. 



imi 

C, 2 _i s 2 

i 

i 

> P 2 R, 2 C t 



FIGURE t. Graphical illustration of two machine problem 


Completion time for job /'on machine I is given by, 

(1) c, 1 - jj sj + p; + r} 

Completion time for job » on machine II is 

(2) C, 2 - (C,L| + 5/ + P,') + P, 2 + R 2 

if C,i | + S, 2 < CJL | + S, 1 + Pf (see Figure 1) 

or 

C , 2 - C,i t + S 2 + P 2 + R 2 
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C, 2 _ I + S, 2 > C,l t + 5,' + P}. 

The combination of (2) and (3) gives 

(4) C 2 = Max [C, 1 - R,' - S 2 . C,l, ] + S 2 + P 2 + R 2 . 

By successive application of (4) using (1), the total elapsed time is given by 

(5) T= C 2 = MaxJ fl Max (5/ - S, 2 + /»,') - JJ (P 2 + R 2 - oj 

+ £ (Si 2 + P 2 + R 2 ). 

The optimal schedule can be found by using Johnson's method to solve a two machine flow- 
shop problem where the processing times of job i on machines I and II are 5,’ — S 2 + P, ] and 
P 2 + R 2 - R,\ respectively. 
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A NOTE ON THE MAXIMIN VALUE OF TWO-PERSON, 
ZERO-SUM GAMES* 


Alvin E. Roth 


Mellon Professor of Economics 
University of Pittsburgh 
Pittsburgh, Pennsylvania 


1. INTRODUCTION 

The theory of two-person, zero-sum games has occupied a special place in the general 
theory of games ever since its introduction by von Neumann and Morgenstem [17], They 
presented the theory of two-person, zero-sum games as a natural extension of the theory of 
rational choice under uncertainty, as modeled by the maximization of an expected utility func¬ 
tion. Their conclusion was that rational players should always play their maximin strategies in 
such games, and should regard the maximin payoff as the "value" ot the game.t 

Subsequent authors have questioned both the validity and generality of these conclusions. 
Ellsberg 14], for instance, points to gaps in the arguments which von Neumann and Morgen- 
stern present in support of their conclusions, and raises questions related to the interpretation 
of the naure of the outcomes of the game, while Aumann and Maschler [1] discuss issues 
which arise in passing from the extensive to the normal form of a game. This latter discussion 
has generated a lively controversy (cf. Aumann and Maschler [2], Davis [3], Owen ]7), Taylor 
1131). McClennen [6] also examines gaps in the arguments presented by von Neumann and 
Morgenstern, and, like Aumann and Maschler, concludes that the prescription that rational 
players should choose maximin strategies cannot be derived directly from the principle that 
rational individuals are utility maximizers. 

There is thus a gap between the rationality assumptions which insure that an individual 
evaluates single person decision problems as a utility-maximizer, and the assumptions which 
insure that he evaluates two-person, zero-sum games in terms of the maximin value. Two dis¬ 
tinct issues are involved here, since the maximin value is intended to be used both to evaluate 
alternative games and to identify maximin strategies as "rational" choices. This paper addresses 
the first of these issues. 

One approach to studying the maximin value is the axiomatic approach of Vilkas (IS) and 
Tijs (14], who consider arbitrary functions defined on games and present axioms which are 


•This work was supported by National Science Foundation Gram SOC75-21820 to the Institute Tor Mathematical Stu¬ 
dies in the Social Sciences, Stanford University, and by National Science Foundation Grant SOC78-09928 to the 
University of Illinois. 

tStrictly speaking, this conclusion applies only to games having saddlepoints in the set of admissable strategies. Such 
games are sometimes referred to as determined games. 
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uniquely satisfied by the maximin value. This paper takes a different, but related, approach, 
and shows that the maximin value of a two-person, zero-sum game can be interpreted as a 
rational individual’s von Neumann-Morgenstern utility for playing the game, and demonstrates 
necessary and sufficient conditions on the individual's preferences over games for this to be the 
case. That is, we will consider conditions which an individuals's preferences must obey in addi¬ 
tion to those which insure that he is a utility maximizer, in order for his utility function to coin¬ 
cide with the maximin value. This approach has been used to study other classes of games, 
e g., cooperative games with side payments (Roth (8), (9]), simple games (Roth [10]), and bar¬ 
gaining games without side payments (Roth [11], [12]). The concluding section of this paper 
will discuss briefly the relationship between the results obtained here and these other results. 

It should perhaps be emphasized at this point that it is not the purpose of this paper to 
further explore the mathematical properties of the maximin value, which are well understood. 
The questions which have been raised in some of the papers cited above concern, rather, the 
interpretation which can be placed on the maximin value as part of a theory of rational (i.e., 
utility-maximizing) behavior. The purpose of this paper is to further explore this latter issue. 

2. THE MAXIMIN VALUE AS A VON NEUMANN-MORGENSTERN UTILITY 

A two-person, zero-sum game g will be denoted by a triple g « (Si,S 2 ,«), where S\ and 
Si are arbitrary sets, and u is a function such that u \S\ x S 2 — R 2 , and such that, for any 
(s.t) € S t x S 2 , u,(s,f) - -u 2 (s,t). The interpretation is that in the game g, the player in 
position / (/ — 1,2) will have available a set of strategies* S„ and will compete for prizes over 
which his preferences are represented by the utility function His opponent’s preferences 
will be precisely opposite his own. (Note that in general it is not necessary that the players 
compete for a single set of physical prizes.+ It is sufficient that the results of the game be such 
that the two players are never in agreement over which of two outcomes is preferable.) 

This formulation defines the game independently of the individuals who will play it. Con¬ 
sequently, it is meaningful to consider the preferences of some individual over the positions in 
a game, or in different games. Letting G be some class of two-person, zero-sum games, and 
letting N = (1,2) be the set of positions, we will be considering an individual's binary prefer¬ 
ence relation P defined on /if x (J, the set of positions in a game. For /, j € N and g , g' € (7, 
the statement U,g)P(J,g') should be read "it is preferable to play position / in game g than to 
play position J in game g'." We will interpret P as a strict preference relation, and define an 
indifference relation I and a weak preference relation W in the usual way.* We will assume that 
the preference relation is defined over the mixture set of lotteries generated by N x G, and 
that it satisfies the rationality conditions necessary for the existence of an expected utility func¬ 
tion v. 


’Typically the strategy set S, could be the set of all probability mixtures of some finite set of pure strategies. 

+Two person, zero-sum games are sometimes interpreted as games in which, for every pair of strategy choices, one 
player’s gains are precisely equal to the other's losses in terms of some physical commodity, such as money. Under 
this interpretation, the game would only be zero-sum in the rare case that the two players have precisely opposed 
preferences for lotteries involving money. In our formulation, the game specifies the utility payoffs to be faced by the 
players. For a given pair of individuals to play a specified game, the monetary rewards to each would have to be ad¬ 
justed in terms of each individual’s utility function for money 

»So, for a.b € N x c. alb if and only if neither aPb nor bPa. and a If b if and only if either aPb or alB. 
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Since we are interested in preferences over a possibly infinite set of alternatives (e.g., let 
G be the set of all matrix games) it will be simplest to refer to the axiomatization of utility due 
to Herstein and Milnor (Si. rather than to von Neumann and Morgenstern’s original axiomati¬ 
zation. Denote by [/>(/,#); (I - p)(j,h)] the lottery which with probability p has an individual 
play position / in game g, and with probability 1 - p play position j in game h. Then we assume 
that the preference relation is defined on a mixture set containing all such lotteries, and that 
this preference relation has the properties of transitivity, continuity, and substitutability 
sufficient for the existence of an expected utility function (Herstein and Milnor [5]). That is, 
there exists a function v (unique up to origin and scale) which preserves preference (v(/,g) > 
v(J,h ) iff ( i,g)P(J-h )) and evaluates the utility of a lottery by its expected utility (v[p(/,g); 
(1 - p)(J.h)] = pv((i,g)) + (1 - p)v(U/i)). 

The fact that we are considering preferences defined over games, which are themselves 
defined in terms of an expected utility function, puts some additional restrictions on the allow¬ 
able preferences. In particular, for any real number k , let g k be the (degenerate) game defined 
by gk - (Si.S 2 ,u) such that |S]| - |S 2 I - 1, and u,(Si.S 2 ) - k. Thus g k is the game* in which 
each player has only one strategy to choose, which results in the player in position 1 receiving a 
utility of k, and the player in position 2 receiving a utility of —k. Since k is an expected utility, 
we know precisely how an individual's preferences should behave over alternative games g k . 
r 'wnally, we assume that the set G contains the set of degenerate games, and impose the fol- 
. /ing conditions on preferences over degenerate games, which have the effect of embedding 
the space of utilities into the space of games. 

(1) (\.gi)PU.g 0 ). 

This simply states that it is preferable to play position 1 in the degenerate game g\ which 
awards the player in position 1 a utility of 1 than to play position 1 in the game g 0 and get a 
utility of 0. 

(2) For any real numbers cand k such that c > 1, 

<!.&)/[- (Uc*); [i - -^j(Uo)] and (i.g 0 )/[j<U*);y 0.*-*)]• 

The first part of this condition states that a player is indifferent between receiving a utility of k 
for certain in the degenerate game g k , or participating in a lottery which gives him a utility of ck 
with probability 1/c and a utility of 0 otherwise. The second part of the condition states that a 
player is indifferent between receiving a utility of k or — fc, each with probability 1/2, or receiv¬ 
ing a utility of 0 for certain. This is precisely what is meant by the statement that the function 
u is an expected utility function. 

Condition (1) permits us to normalizet the utility function v in the natural way, and we 
henceforth take v(l„g 0 ) - 0, and v(l,^ f ) - 1. Condition (2) then completes the embedding of 
utilities into the space of games, giving us the following result. 

PROPOSITION 1: For any real number k, v(l,g*) — k. 


•Technically speaking, our definition allows there to be more than one game g k , since we could allow different one- 
element strategy sets. This difference is inessential, and can be ignored. 

(Recall that v is determined only up to choice of origin and scale. 
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PROOF: Since v is a utility function, the first part of (2) implies that, for c > 1, 
- vll/cO.fc.*); (1 - (l/c»<Uo)]. But v is an expected utility function, so the utility 
of a lottery is its expected utility; ie., v[l/e(l,&*); 0 - (l/c))(l,go)] “ (l/c)v(l,g rt ) + 
(1 - (l/c))v(l,g 0 ) ” 0/e)v(l,&*). Thus, for any k , v(l,g rt ) - cv(l,g t ) for c > l. For 
0 < c < 1, let d - 1/c, and k' - ck. Then vO,®*) - dv(\,g k ), so v(l,g ft ) - cv(l.g*) for 
any k and any c > 0. In particular, for any k > 0, v(l,g*) - fcv(l,gi) - k. But the second 
part of (2) implies that v(l,g*) - -v(l,g_ t ) for any k , and so v(l,g*) - k for all k , as was to 
be proved. 

So far we have imposed conditions on preferences over degenerate games which reflect 
the fact that they are associated with the underlying expected utility function, but we have yet 
to impose any conditions which reflect the zero-sum nature of all games in G. The following 
conditions will suffice for our purposes. 

(3) For all g € G, <U 0 >/[y 0.#); y (2,*) 

This condition states that a player is indifferent between getting a utility of zero for certain in 
the degenerate game g 0 or participating in a lottery which gives him an equal probability of 
playing either position in any two-person, zero-sum game g. The motivation for this require¬ 
ment is the idea that if a given zero-sum game g yields an advantage to a player in one of the 
positions, then it must yield a corresponding disadvantage to the player in the other position. 
An immediate consequence of condition (3) is the following: 

PROPOSITION 2: For all g € G. v(l,g) - -v(2,g). 

PROOF. (3) implies that v(l,g 0 ) - v[|/2(l,g); 1/2(2,g)] - (l/2)v(l,g) + (l/2)v(2,g). 
But v(1 ,g 0 ) - 0, so v(l,g) - —v(2,g). 

To state the next condition, consider a game g « ( Si,S 2 ,u) and a given position i € N. 
for every strategy s € S„ define k(s) - k(g.i.s) - inf u,(j,r), where j * /.* Then in addition 

f€S / 

to the conditions already imposed on the preferences, we require the following. 

(4) For any g € G. i € N. and s € S„ (i,g) WUgki,)). 

If the game g were such that 5, contained only a single strategy, then condition (4) would sim¬ 
ply be a version of the "sure thing” principle, which states that any prospect is at least as desir¬ 
able as the worst outcome which can result from that prospect. As it stands, the condition also 
reflects that individual i is free to choose any strategy s in 5,. Condition (4) implies the follow¬ 
ing proposition. 

PROPOSITION 3: For any g € G, 
v(U) ^ su|i inf Ui(s.r), 

and 

v(2,g) > sujj inf u 2 (/.s). 


*So t(g. I.s) • inf U| (s.t), and k(g,2.s) m inf uj(r.s). 
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That is, the utility of playing a game is at least as large as the "security level" which the game 
provides. No further assumptions are needed to permit us to reach the following conclusions. 

PROPOSITION 4. For any g € G, 

inf suj> W|(s,/) > v(l,#) > suj» inf u,(s,r). 


PROOF: The right hand inequality is simply the first part of Proposition 3, while the left 
hand inequality follows from Proposition 2 and the second part of Proposition 3, together with 
the fact that uAs.t) - -u 2 (s,t). 

A game g is said to have a saddlepoint if there exist strategies s 6 S, and t € S 2 which 
achieve the inf sup and sup inf of Proposition 4, and make them equal. For such games. Pro¬ 
position 4 has the following immediate consequence. 

PROPOSITION 5: If* € C has a saddlepoint, then 

v(l.j) - max min uAs.i) — min max ui(s,r). 
s€S, HS } r€Sj j«S, 

The famous Minimax Theorem of von Neumann [16] establishes that if g is the mixed exten¬ 
sion of a finite two-person, zero-sum matrix game, then g has a saddlepoint and Proposition 5 
applies. 

DISCUSSION 

The previous section establishes conditions on an individual's preferences over games 
which yield the maximin value as his utility for a given game. Insofar as these conditions are 
consistent with our other notions about zero-sum games, this shows that the maximin value is 
consistent, both formally and substantively,* with the notion of rationality embodied in utility 
maximization. The "machinery" necessary to establish this result shoiuld not, however, be per¬ 
mitted to obscure the essential simplicity of the argument. 

The two substantive conditions on preferences are (3) and (4). Condition (4) serves to 
establish a lower bound on the utility of playing the game from either position. Condition (3), 
by fixing the utility for one position in the game with respect to the utility for the other posi¬ 
tion, combines with (4) to produce an upper bound on the utility. For a game with a 
saddlepoint these bounds coincide, and so the utility is determined. 

Since the maximin value is often interpreted as being an unnecessarily conservative 
assessment of a game (cf. Etlsberg's [4] "reluctant duelist"), it is interesting to note that condi¬ 
tions (3) and (4) actually prevent the utility function from reflecting an excessively conserva¬ 
tive attitude. In particular, if g is a game without a saddlepoint, then the utility of playing g in at 
least one of the positions must be strictly greater than the maximin value for that position. 
This is because Proposition 2 requires that v(I,^) — — v(2,#), and, in a game without a 
saddlepoint, the maximin values for each position are not the negatives of one another. 


‘Any numerical index for certain events can of course be extended to a function which preserves expected value over 
lotteries. The purpose of this paper is to investigate what substantive assumptions about preferences need to be made 
in order to interpret this extension as an expected utility function. 
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Note that if the class G does include games which have no saddlepoint, then the condi¬ 
tions imposed in the previous section do not completely determine the utility function v. To 
do so, it would be necessary to precisely specify a player's attitude towards the additional uncer¬ 
tainty present in games with no saddlepoint. The same phenomena occurs in studying coopera¬ 
tive games, in which it is necessary to specify a player's attitude towards strategic risk (cf. Roth 
191). 


As a final observation, note that, although the maximin value is linear in degenerate 
games, it is not linear in arbitrary games. Let g and h be finite matrix games of the same 
dimensions, let p be a probability (i.e., p € 10,11), and let g' - pg + (1 - p)h (where we 
employ the usual conventions of matrix arithmetic). Then each element of the matrix g' is the 
expected value of the corresponding element of the random matrix defined by the lottery 
[/>*;( 1 - p)h\. We might conjecture that a player would be indifferent between playing a given 
position in the game g' or taking the same position in the lottery. But for preferences obeying 
the conditions considered in this paper, this is not the case. That is, v[p(i,g); (1 - p)(i,h)\ 
equals pv(i.g) + (1 - p)v(i,h) and is in general not equal to v(i.g'), for / - 1,2.* In this 
respect, the preferences considered here over zero-sum games resemble the preferences over 
bargaining games considered in Roth [111, (121 more than the preferences over games with 
sidepayments considered in Roth [81, [9]. 
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ABSTRACT 1 

Certain types of communication nodes can be viewed as multichannel \ 

queueing systems with two types of arrival streams. Data arrivals are character¬ 
ized by high arrival and service rates and have the ability to queue if all service 
channels are busy. Voice arrivals have small arrival and service rates and do 

not have the ability to wait when the channels are full. Computational pro- I 

cedures are presented for obtaining the invariant probabilities associated with 
the queueing model 


1. INTRODUCTION 

The steady growth of the data processing industry and the telephone network over the last 
four decades has introduced computer-communications networks as efficient transportation 
vehicles for the remote sharing of information data bases. Recent Defense Communications 
Agency studies have shown the desirability of a network which integrates voice, interactive, and 
bulk data for the 1980's (Rosner [10], Schmitz, Saxton, Huang and White [11]). Several addi¬ 
tional studies relating to information processing growth in the next few decades portend new 
data/voice services with substantially increased data flows (Rich and Schwartz [8], Rosner [9]). 
There are several different time division integration methodologies proposed for combining 
voice and data demands over the available channel bandwidth (Coviello and Vena [2], Fischer 
and Harris [4], Frank and Gitman [5]). A competitive allocation scheme allows the data and 
voice calls to compete for time slots using a first-come, first-served scheme. 

Simulation models have been developed to study competitive allocation schemes; how¬ 
ever, because of their expense and the fact that rare events are of interest, simulation studies 
are not necessarily appropriate. Analytical models have also been developed to describe com¬ 
petitive allocation (Bhat and Fischer [1], Fischer [3], Fischer and Harris [4]). The analytical 
models developed to date all share the common drawback of not being tractable for problems of 
reasonable size. The purpose of this paper is to present an analytical model of the dynamic 
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movement of voice and data calls within an integrating device and give a tractable solution pro¬ 
cedure for the model. The allocation scheme is that used by Bhat and Fischer [lj. The 
methodology for the solution scheme combines procedures developed by Neuts [6,71 and 
Wong, Griffin, and Disney [121. 

2. THE MODEL 

There are two classes of calls that arrive to the service center: data calls and voice calls. 
The arrival process of data calls is assumed to be Poisson with rate A.) and the arrival process of 
the voice calls is assumed Poisson with rate X 2 ; the two arrival processes being independent. 
Although each data call actually consists of a random number of discrete bits, the length of 
each data call is assumed to be adequately approximated by a continuous exponential random 
variable with mean l/p, and the length of the voice call is assumed exponential with mean 
\/n 2 The characteristics of data calls are frequent arrivals and short length whereas the voice 
calls have infrequent arrivals with a long service time. A buffer is available that can store 
incoming data calls when all channels are busy. The size of the buffer is such that no effective 
limit is placed on the queue size. Voice calls, however, may not form a queue so that any 
incoming voice call is lost whenever all channels are busy. The service center has c channels 
using a FIFO discipline. 

The data/voice integrating system is modeled as a Markov process with a two dimensional 
state space given by 

£ — {(/>,/»): m =0,1.c and n =0,1, ...}. 

The vector p (ordered lexicographically) denotes the steady state probabilities where p(n,m) is 
the steady state probability that there are n data calls in the system and m voice calls in the sys¬ 
tem. For ease of notation the vector p is partitioned as 

(1) P - (Po- Pi, ■•• ) 

where the mth component of p„ is p(n.m). The vector p is the solution to the equation 
P(? - 0 

(2) Pl=l 

where 1 and 0 are vectors of all ones and zeros, respectively, and Q is the transition rate matrix 
for the Markov process representing the data/voice system. The matrix Q is infinite and is 
given by 


(3) 


Q 



O 


M c A c A 
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where the submatrices are defined, for i.j.n -0.1.c, by 

( X. if» - j. 

(4 > -10 otherwise: 

fl/i i if i - j and j < c ~ n. 

(5) M„U,j) - (c - j)m if i - j and j > c ~ n, 

0 otherwise: 

k 2 if / - J - 1 and i < c - n - 1, 

M2 if i - j + 1, 

a„(i) if i - j. 

0 otherwise; 

where a n (i) is such that the row sum is zero. 

3. SOLUTION PROCEDURE 

The general structure of the matrix Q is identical to the type investigated by Neuts [6) 
and similar to the finite system investigated by Wong, Griffin and Disney 1121. The solution 
method as given by Neuts [61 is a two step process. First, it is necessary to find the matrix R 
such that 

(7) R 2 M c + R A l + A - 0. 

The matrix R gives the relationship 

(8) p f+ * - p c -\R k+ ' for*-0,1. 

The second step is to let p = (p 0 .p f _i) and then solve for p u'i'-g equations 

(9) p T - 0 
and 

Pl + P ciRU-R)-' I- 1 

where 


(6) A„(iJ) " 



It will turn out that for this communication problem, the matrix R of equation (7) is easy to 
obtain whereas 0 from equation (9) gives computational difficulties which will be overcome by 
using a technique utilized by Wong, Griffin and Disney [12]. 
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The matrix R is lower triangular since the matrices M e , A c , and A are lower triangular. 
Therefore, the diagonal elements of R 2 are simply the square of the diagonal elements of R. 
Thus, the following equation results when considering the diagonal elements in (7) for 


/ = 0.c 

(10) r,} (c — /)/*| — r„Ui + ifi 2 + (c — /)/i|) + A| — 0. 
The solution of (10) yields 

r n — A |/(cm 2 + Aj) 
and for i - 0. c - 1 

r„ — U| + ifij + (c — i)fi | — l(X i + ifi 2 + (c — / )/* i) 2 

(11) —4a,(c - /)mi1 i/2 >/(2(c - 


The off-diagnonal elements are not quite as straight forward but by considering the /,./ element 
with / > j, equation (7) yields 

(C - ,/')m| £ r ik r ki + £ r 'k A t (kJ) - 0 
which in turn yields 

(c - j)fi\ [r,>(r„ + r n ) + j) r ik r kl ] 

+ 0 + 1 r r./+t ” U, + jfii + (c — j)n\\rij 
and thus for / > j 

(12) r ti “ (O' + 1 >M 2 r u+\ + (c —j)n i X + jfk 2 + (c—j)fi | (1 - r (/ — r (/ )) 

where the sum in equation (12) is defined as zero if j - / - 1. Equation (11) and (12) give an 
easy iterative procedure where the diagonal elements of R are first computed, then the ele¬ 
ments such that j - i — 1 are computed then the elements such that j - / - 2, etc. In this 
manner the exact expression for R is obtained and the solution of equation (9) remains. 

The obvious difficulty with (9) is that for a typical system in which c - 48, the dimension 
of fis 2352 x 2352. In order to obtain a solution to (9) it will be reduced to solving a 49 x 49 
system. 

Equation (9) can be rewritten as 

(13) Po^o + Pi A/1 - 0 

(14) P*-|A + Pn-4* + p* + iA/* + i - 0 for k - 1.c - 2 

(15) p f . 2 A+p c _,(/4 c _, + «M ( )-0. 

Each of the submatrices in the above equation are of dimension (c + 1) x (c + 1). New 
matrices B k are defined with dimension (2c + 2) x (2c + 2) by 

-A k A' 1 / 1 

-M k+] A' 1 0] for * “ 1 . C “ 2 ' 


*- 
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II now follows from (14) that 

(16) <Pi-i.P*>- <P*.P* + i)B* 

and thus 

(H) (Po.P,) - (P.-2 .Plli- 

Combining equations (13) and (15) with (17) results in 

[ ^ol 

(18) p r _,(U r _, + RM r ]\-'.l)B ^ 2 ... fl, ^ I - 0. 

Equation (18) gives a solution for p, . t that is unique up to a multiplicative constant. Any such 

solution is obtained and equation (16) is used to obtain p ( ._ 2 , p r _j.Po The norming 

equation of (9) is then used and the steady state probabilities are determined. 

A unique combination of two previously determined solutions techniques thus yields a 
computational procedure that can be used to solve a class of problems in the telecommunica¬ 
tions field. Typical measures of effectiveness can easily be determined in the same manner as 
in Neuts [6). Thus, the utilization of analytical models for such telecommunication systems is 
now feasible. 
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ABSTRACT 

This paper attempts to resolve the existing confusion concerning missing 
operations. Scheduling problems are classified in two groups: (i) null- 
continuous (NO—comprising the problems where an optimal schedule remains 
optimal on replacement of arbitrarily small processing limes (existing opera¬ 
tions) with zeros (missing operations): (ii) null-discontinuous (NDC) — 
comprising those problems which are not null-continuous. 


A "zero processing time” of an operation refers to either of the two following contingen¬ 
cies: (i) An actual operation whose processing time tends to zero. Thus, if t„ denotes the pro¬ 
cessing time of operation 0, in job J„ then for a sufficiently small positive number e, scheduling 
problems with t 0 - e and fy — 0 have the same optimal schedules (algorithms), and e may rea¬ 
sonably be replaced with zero to facilitate calculations, and (ii) A nonexisting (missing) opera¬ 
tion. 


Since an infinitesimal (arbitrarily small) processing time operation (i) has a starting time 
while a missing operation (ii) has not, the two types have a different effect on a scheduling 
problem and must be differentiated to prevent ambiguity. Accordingly, we propose to designate 
an infinitesimal processing time operation as e, and a missing operation as a zero. 

Scheduling problems involving operations with arbitrarily small processing times (existing 
operations where, for all practical purposes, the length may be considered as zero, but have 
starting times) are basically the same as those with the usual strictly-positive processing times 
and do not merit separate consideration. 

Whether a discipline (flowshop or openshop) allows missing operations or not depends on 
its definition, without clear preference for one definition over another. However, for the sake 
of understanding, the definition (whatever it is) must be known and accepted. We propose to 
define flowshop and openshop disciplines allowing missing operations. (Note, [2] and [3] 
define flowshop allowing missing operations while 14], [9] and [14] implicitly assume that 
flowshop does not allow missing operations.) Accordingly, For 0 in the notation n/m/y/b*, 
y € |F,0) comes under that definition and the processing time of the 0 1 operation of job 
7,, (i • 1,2, ... , w, 7 - 1,2, .... m), is a nonnegative integer, t 0 ^ 0—positive if the opera¬ 
tion exists and zero if it does not. In the case of a flowshop. or an openshop, where all njobs 
have m operations—i.e., where the processing times of all mrt operations are strictly positive— 

*The notation is that proposed by Lenstra [91 and Rinnooy Kan [14], 
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the constraint t u > 0 V ( i,j ) is added, and the designation becomes n/m/y. t„ > 0/8, 
v 6 (F.0). 

Panwalker, Smith and Woollam [Ul, demonstrated that the commonly held view that 
there exists an optimal permutation schedule for n/3/F/C ma% and n/m/F.\ no-wait/C ma ,, 
m > 3, turns out to be correct only for positive processing times.* (The latter fact was recog¬ 
nized by Baker (21.) 

1. NULL-CONTINUOUS AND DISCONTINUOUS SCHEDULING PROBLEMS 

The general rule is that optimal schedule remain unchanged under small changes in the 
close neighborhood of the given processing times, providing there is no replacement of t with 
zero (deletion of an infinitesimal operation). Where such replacement is resorted to, some 
problems remain unaffected while in others the optimal schedule undergoes total transforma¬ 
tions; in the latter category, an optimal algorithm and schedule for > 0 V(/,y) are not 
optimal for t u > 0. Moreover, there are problems for which there exists an efficient optimal 
algorithm (belonging to P) for the first case (/ 0 > 0 V (#,y)) while for the second (t„ > 0), the 
problem is NP-Complete. 

DEFINITIONS 

Null-continuous (NC) scheduling problem — one where an optimal schedule remains 
optimal on replacement of arbitrarily small processing times with zeros and vice versa. Thus, 
an optimal algorithm and schedule for t 0 > 0 V (ij) are also optimal for t„ > 0. 

Null-discontinuous (NDC) scheduling problem — one which is not null-continuous. 
Thus, an NDC problem defines two different problems, one with strictly positive processing 
times (/y > 0 V (/j)—missing operations are not allowed) and the other where zero processing 
times are allowed (f /y > 0—missing operations are allowed). 

Informally, a problem where all operations with arbitrarily small processing times for 
every machine can be shifted to the beginning or the end of the schedule without affecting the 
measure of performance or violating the constraints—is NC. To show that a problem is NDC, 
it suffices to work out an example where replacement of arbitrarily small processing time with a 
zero results in a different optimal schedule. 

The following numerical examples of NC and NDC problems demonstrate the fact that an 
NDC problem defines two different problems—one where > 0 and the other where t 0 ^ 0. 

EXAMPLE 1: 

NC instance — 5/2/F/C ma J with processing times as per Table 1. 


•[Ill was brought to our attention by the referee. At the time of writing the first version of the paper, we were 
unaware of 111] and independently reached the same conclusions. 

+The n/2/F/C mn and n/2! F, no-wait/C m „ problems are NC and NDC, respectively. 
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TABLE 1 — Processing Times 


Job 

OperatiofN, 

(machine) 

1 2 3 4 5 

1 

2 

4 e e 3 2 

2 2 2 « 4 


The problem is NC since arbitrarily small processing time operations on the first and second 
machines can be shifted to the beginning and the end of the schedule, respectively, without 
violating any constraint or affecting C max . Thus, optimal schedules constructed by Johnson’s 
algorithm [8] remain optimal after replacement of « with zeros. 

NDC Instance—5/2//", no-wait/C max with processing times again as per Table 1. For 
tjj > 0 V (i,j) Gilmore and Gomory's algorithm [5], [13] yields an optimal schedule as 
presented in Figure 1. The algorithm itself is not affected by replacing e with zeros, but the 
resulting schedule is no longer optimal (Figure 2). 



nr 

— m 



J 2 J 3 

J 5 

M 1 

0 2 l 

6 

10 11 



NP-Completeness of a scheduling problem is proved by reducing a known NP-Complete 
problem to an instance of the problem in question that in many cases contains zero processing 
times. In the light of the preceding discussion only if the problem in question is NC then 
results obtained for t 0 > 0 V(/J) are also true for ^ > 0, and vice versa; if NDC, the problem 
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may refer e'ther to t„ > 0 V(/,./) or /„ > 0 with a different approach and specific optimal algo¬ 
rithm and schedule in each case. Thus, before zero processing times are allowed, the problem 
in question must be characterized as NC or NDC. 

2. CLASSIFICATION AND COMPLEXITY OF NC AND NDC PROBLEMS 

A first attempt has been made to classify some of the scheduling problems into NC and 
NDC, and to examine the complexity of each problem (the results are obtained in Table 2). 
As explained previously, each NDC problem defines two different problems, one with strictly 
positive processing times (/„ > 0 V(/,./')) and the other with nonnegative ones, (/„ > 0). The 
fact that for a particular problem, the case where /„ > 0 V(/,y) is a special case of r„ > 0 and. 
thus, its complexity is lower than or equal to that of the latter (theoretically speaking, all com¬ 
binations are possible) is demonstrated in Table 2. 


TABLE 2 — Classification and Complexity of NC and NDC Problems 



Problem 

Complexity 

t„ > 0 

^ 0 


n/2/F. no-wait/C max * 

0(/t 2 ). [5] 

Unarv NP-Complete, [151 


n/3/K no-wail/C max 

?t 

Unary NP-Complete, [4] 


n/4/F, no-wait/C max 

Unary NP-Complete, [12] 

Unary NP-Complete. [12] 


n/2/F, r, > 0/C max 

Binary NP-Complete, [10] 

Binary NP-Complete. [10] 


n/2/F tree/ C max * 

Binary NP-Complete. [10] 

Binary NP-Complete. [10] 


«/2//vtree'/C max § 

0(n log n), 116} 

? 


n/3/ Ff C max 

Unary NP-Complete, [10] 

Unary NP-Complete. [4] 

u 

n/2/0, no-wait/C max * 

Unary NP-Complete, [15] 

Unary NP-Complete. [15] 

z 

n/2/J, n, - 2, no-wait/C max * 

Unary NP-Complete, [15] 

Unary NP-Complete. [15] 


n/2/F/ZC, 

Unary NP-Complete, [4] 

Unary NP-Complete. [4] 


n/2/F no-wait/I C, 

?+ 

Unary NP-Complete, [4] 


n/2/0, no-wait/I C, 

? 

Unary NP-Complete, [1] 


n/ 2/0/EC, 

? 

Unary NP-Complete, [1] 


n/2/J/nj -= 2, no-wait/LC, 

? 

Unary NP-Complete, [4] 


/7/l/seq.dep./C max 

Unary NP-Complete, [14] 

Unary NP-Complete.[14] 


(Travelling salesman) 




n/2/F/C m , x 

0(« log n), [8] 


u 

n/2/J, n, < 2/C max 

0(r» log /?), [7] 


z 

n/m/0/h 

Depends on the number of machines (m) and the regular] 

_1 


measure of performance (8). ; 


•Flowshop. jobshop and openshop disciplines are defined allowing zero processing times. 

+ Prize carrying open problem. |6|. (9|, |I0|. 

tlree—Tree precedence relations where job J, may start only after job 7, has been completed. 

§tree' - Tree pr cedence relations where, on each machine, job J, may start only after J on this machine 
has been completed. 
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